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MATHEMATICAL FORMULATION OF CARDIOVASCULAR 
DYNAMICS BY USE OF THE LAPLACE TRANSFORM 


SIDNEY ROSTON 
DEPARTMENT OF MEDICINE, SCHOOL OF MEDICINE 
UNIVERSITY OF LOUISVILLE 


By means of the Laplace transform, the behavior of a simplified model 
of the cardiovascular system is mathematically formulated. This formu- 
lation allows mathematical expression of the periodicity of the cardiac 
output and the systemic response. With the cardiac output represented as 
half of a sine function cycle, the systolic aortic pressure becomes the 
sum of a sine term and exponential terms, while the sum of the expo- 
nential terms alone represents the diastolic pressure. The character- 
istics of the mathematical expressions for systole and diastole are 
analyzed, and some relationships of potentially practical value are de- 
rived. Variation in the parameters of the system yields mathematical 
results consistent with the expected physical ones. 


Introduction. The use of the air chamber (Frank, 1899) as a hydro- 
dynamic model for the mathematical formulation of cardiovascular dy- 
namics has received a great deal of attention over the years. Despite 
a number of modifications, this model has not yielded satisfactory 
mathematical results (Aperia, 1940). In retrospect, it now appears 
that the inadequacy may, in part at least, have been due to the 
weaknesses of the mathematical techniques rather than to the in- 
sufficiency of the model itself. No other completely satisfactory 
mathematical approach to cardiovascular dynamics has been formu- 
lated. Recent contributions to this formulation have included two 
basically pragmatic techniques, namely, the ‘‘pressure-pulse’’ con- 
cept (Hamilton and Remington, 1947) and the use of regression 
equations based upon post-mortem simulation of systole (Starr e¢ 
al., 1954). Each of these approaches has been unfavorably criti- 
cized (Brotmacher, 1957; Bridwell e¢ al., 1956). 

The simplicity of the chamber model makes desirable its recon- 
sideration with a mathematical technique not heretofore used, 
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namely, the Laplace transform (Gardner and Barnes, 1942). The 
Laplace transform has had widespread application in engineering 
calculations but has not to any significant extent been applied to 
physiological problems. The technique is particularly useful in 
situations involving periodic functions of unusual shape, such as 
the cardiac output. The mathematical solution of a specific prob- 
lem in cardiovascular dynamics presented in this work comprises 
also a demonstration of the potential usefulness of the Laplace 
transform approach in a wide range of physiological problems. 


Theoretical results. Certain assumptions facilitate the mathema- 
tical formulation of cardiovascular dynamics. A criterion for the 
reasonableness of these assumptions is the usefulness of the 
theoretical results. The cardiac output W(¢) is, for mathematical 
simplicity, approximated as a pressure-independent function of time 
’ ¢ by the positive half of a sine function cycle, the flow extending 
throughout systole (Figures 1 and 2), The isometric contraction, 
as well as relaxation, of the ventricles is included in diastole 


FIGURE 1. Diagram of the chamber. 


rather than systole. Blood is assumed to be a homogeneous non- 
compressible fluid. The peripheral circulation is represented by a 
single effective resistance FR, the value of which is not a function 
of time but rather of such factors as the degree of arteriolar con- 
striction, A single chamber with elastic walls represents the 
aorta. Whether or not the volume of the theoretical chamber is the 
same as, or indeed bears any fixed relation to, the actual volume 
of the aorta remains to be determined by practical application of 
the theory. Pressure changes are assumed to be instantaneously 
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FIGURE 2. Representation of the cardiac output W(¢), with m= 4, a= 


1 


mi z 
0.8 seconds, pee 15.7 sec. The solid lines represent W(t). The 


broken lines represent components of W(¢), as follows: Curve 1: 
bu(t) sin 15.7%; Curve 2: bu(¢—0.2) sin 15.7(¢— 0.2); Curve 3: bu(t— 
0.8) sin 15.7(¢ — 0.8); Curve 4: bu(é— 1.0) sin 15.7(¢ — 1.0). 


transmitted throughout the chamber, gravity and standing waves of 
pressure in the distal aorta being neglected (Hamilton and Dow, 
1939). The pressure in the aorta near the aortic valve is taken as 
the chamber pressure, P(t). The total volume of the chamber is 
the sum of the constant volume V, at which pressure first appears 
and the volume V (¢) in excess of V,. With K the constant of elas- 
ticity of the aorta, it is assumed that the relation between P (¢) and 
V(t) has the form, V(¢) = =, It is interesting to note that a 
linear relation between pressure and volume exists in the aortas of 
cadavers within a fairly wide range of age and pressure (Reming- 
ton, 1948). The relation between the flow through the outlet, F'(¢), 
and the pressure is postulated to have the form F(t) = P(t)/R. Al- 
though this linear relation does not strictly hold true for the flow 
through a capillary bed (Green et al., 1944), it may be reasonable 
for the body as a whole. The effects of closure of the aortic valve, 
such as the dicrotic notch in the pressure tracing and the brief flow 
of blood in the aorta toward the heart (Green e¢ al., 1940), are neg- 
lected. In the steady state, the average output from the chamber 
equals the average input, and there is no net accumulation of fluid. 
The effect of the venous pressure in the body is neglected. 

With these simplifications, the basic equations assume the form 


dV (t)/dt =W(t)-F(t), (1) 
ca + 2 =W(t), since F(¢) = =” and - 


dV(t) 1 dP(t) 
eee K eit 
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Then, in terms of the Laplace transform (Gardner, Loe. Ctt.), 


P 
Es [sP(s) - P(0)] + 2 =W(s), where P(0) = value of P(¢) at 
K 
t = 0, and 
KW (s) P (0) 
= —————_— + ——_.. 3 
AO! oo K/RO SS ERPR (3) 
Now, the equation for W (¢) in Figure 2 is 
n mn 
= bul(t — ine (i 
W (t) Zs | u(t — na) sin (t — na) + 
(4) 


a _ mn a 
bu(e—na~<) sin RE (na 2)), 
m a m 


where u(t¢-—na)=0 for t<na; 
=1 for t2na; 


a a 
u{t-na-—}]=0 for t<na+—-; 
m m 


=f for ¢2 na + 23 n being a positive 
integer or zero. 
Here, a is the period and 1/a is the frequency of the cardiac cycle, 
with flow from the left ventricle into the aorta being actually pres- 
ent only during systole. Systole is the initial 1/mth portion of each 
cycle, m an integer. The function representative of W (¢) during the 
first cycle, W, (2), then is 


W(t) = (sin na) u(t) + bu ( _ | sin a ( <i =), (5) 


The Laplace transform of 


. mrt. bam 
bsin —— jis ————- 


a meets 
a\s? 4+ —— 
a 


The Laplace transform of W, (¢) then is 


br - 4s 
W,(s) = a (, +e ™ ), (6) 
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Now, the inverse transform of the term 
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_, mak ; 
where wy = tan ies The inverse transform of 
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Then, 


_ Ke K bam : - Rena) 
P(t)=P(0)e *+ ep Ee —na)e + 
a\ — + n=0 

(5 a ) 


_x t—na—— 
u(e-na - = Je RI |. 


m 
4 _ (9) 
Wane |) juce ~— na) sin S (t — na) - uf + 
tae \ a0 
( : <) ae 


st-ne- 2) anf e-me- 2) 


; = . mart 
The cardiac output per cycle, D, is equal to | B sin reo dt= 
) 
2ab 
=. Len, 
Mt 


7 Dam 
Qa 


(10) 


Since the volume of the chamber is the same at the end of each 


[ Pwd. But 


cycle in the steady state, D -{ F(t) dt = 
) ) 


av 


[ P(t) dt=aP,,,, where P|, = average pressure. Hence, 
0 


D waren ee (11) 


Discussion. The expression for P(¢) in (9) contains essentially 


three components, There is, in the first place, a non-repetitive term 
_Kt 


P(O)e R which represents the transient response and gradually 
fades to zero. In the second place, when ¢ = na and ¢ = na + a/m, 
-  (t-na) - § (t-na-*) 

and Fe m’, respec- 
Kbnm 


“ =) 
ra eae 
R? a? 


exponential terms equal to Fe 


tively, are added, where F = (Figure 3). Finally, 
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FIGURE 3. The curve of chamber pressure for m = 4, a = 0.8 seconds, 
Afi. = second !, K = 1600 dynes/cem5, D=72 ml. Curve 1: A sine 


term, in mm Hg: 43.2 sin $15.7(¢t — 0.8) — 1.51. Curve 2: An exponential 


term, in mm He: 43.2 e—*. Curve 3: P(t), in mm Hg, with omission of the 


term P (0) e R , The broken lines represent the sum of exponential terms 
during systole. 


at times ¢ = na and t = na+a/m, sinusoidal terms with amplitude 


Kb 
________ and phase angle - w~ appear (Figure 8). Each 


K2 mn? > 
tee ae) 
sinusoidal term starting at ¢ = na produces an effect for a time a/m 
and is cancelled thereafter by the sine term appearing at ¢= na + 
a/m. In the steady state, the systolic phase of the volume and 
pressure is the sum of a sine term and exponential terms, while the 
diastolic phase is the sum of exponential terms alone. 

The exact mathematical handling of the systolic phase is diffi- 
cult, partly because of the presence of the phase angle - yw. The 
peak of systolic pressure would be expected to occur at a different 
time from the peak of the sine function because of the presence of 
the decreasing exponential terms during systole. The pulse pres- 
sure PP is nearly equal to twice the amplitude of the sine term 
less the decay P, in the sum of the exponential terms during the 7 
seconds between the onset of systole and the peak of systolic 
pressure (Figure 3). The decay in the exponential terms is de- 
termined up to the time of the peak of systolic pressure rather than 
the maximum value of the sine term because the pulse pressure is 
measured in terms of the systolic pressure peak, P, can be closely 
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approximated by the drop in pressure during the first ¢ seconds of 
diastole. The approximation will be exact if m = 2, since, then, 
the sum of the exponential terms will have the same value at the 
beginning of systole and diastole. If m is greater than 2, the sum 
of the exponential terms at the onset of systole will be less than 
that at the onset of diastole, since the decay during the shorter 
systole will be less than the decay during the longer diastole. 
The approximation will, even in these circumstances, be a close 
2 se 
one in practice, however. Then, if pa is much less than os ; 


which is ordinarily the case in the cardiovascular system, the 


KD 
amplitude of the sine terms becomes aL. and, as a result, ae 


(PP + P,)/2 and 
ed ae we 
D.2 


RK (12) 
The steady-state diastolic pressure may be expressed as 
SAKE n Kna 
P, (t) = Fe Lape Kc -naje ® 4 
ee (13) 


u(t = Te =) | . 
m 


The pressure during the first diastole has the value 


— Kt a\ Ke 
Fes [vos u(e~ 2) eh], 
m 


that during the second diastole, the value 


- Kt a\ Xa Ka 
Fe uy ru(e- eR sue ayer + 
m 


U (! -a- a) risa) 
m 

aw Ke 
and so on, tha transient term P(0)e * being neglected. The 
mathematical expression for each diastolic pressure consists of 


components independent of time and constant during a given dias- 
— Ke 


tole, and one term, e *®, which is a function of time. Strictly 
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speaking, no diastolic pressure pattern is exactly like any other, 
since, theoretically, the exponential terms added in each cycle 


only reach zero value at ¢ equal to infinity. Eventually, however, 


he nKa KR (ae2) 
e becomes negligible in comparison with e ®; e®* ™ 


eet © Pen eee) 
negligible in comparison with e ™*; and so on, and a repetitive 


pattern is reached. Since the volume and pressure at the end of 
each diastole are, in practice, considerably more than zero, the 

K a 
=~ ne) 


becomes 


term e must retain a significant positive value at least 

at the time ¢=na+a. For each diastole, then, the pressure is ex- 
Kt 

pressible in the form P, (t) = Cpe R where Cy is a constant 

within each diastole, but increases with each additional diastole 

until a stable value is reached, and ¢ extends: continuously from 

zero to infinity. Then, 


Kt 
InP=InCy-—, (14) 


and a graph of In P as a function of ¢ is a straight line with the 
slope — K/R(= — G@) for each diastole, with the changes in C’ hav- 
ing no effect upon the slope. A combination of equations (12) and 


K 
(11) with Bia G yields 


(15) 


adP(t 
It is interesting to note that if b= 0 and W(¢) = 0 in (2), = ) = 


K 
~ a? and In P is again a linear function of ¢. Only by the use 


of the Laplace transform, however, has a mathematical representa- 
tion of the diastolic behavior of the model in the steady state been 
obtained. The basic similarity of this steady-state result to the re- 
sult when W (Z) is zero is expected and indicates the validity of the 
Laplacian representation of diastole. The steady-state arterial 
blood flow in a limb considered as an isolated segment of the body 
cannot be represented mathematically in the same way as that in 
the aorta. This results from the fact that the function W (¢) for the 
limb, being part of the output from the aorta, cannot reasonably be 
approximated by the value zero during diastole. 
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Wherever K and R appear in the expressions for V (¢), they occur 
in the relationship K/R. Thus, a change in K will be balanced by 
a proportional change in F# as far as the volume is concerned. How- 
ever, even in the event that & changes proportionately to K, the 
pressure will change in the same direction as K, since P(¢) = 
KV(t). Variation of V(¢) and P(t) with changes in the parameters 
is of interest. As R becomes infinitely large, V (¢) approaches the 
limit 


ViO\n ss = — + =~)" juCe ~ an) ru(e~an == aN 
n=0 
u(t + na) sin {ete ~na) - a (16) 


a 2 7m a 7 
u(t —-na — — }) sin y—(t-—na - — - =} : 
m a m 2 
The second of the three components in this expression contains a 


series of unit step functions multiplied by one-half of the stroke 
output, a new step function appearing each time ¢ becomes equal 


to na and na+a/m. In addition, ~ = tan ! « = = and the sine 


term reaches its maximum value at the end of each systole. Thus, 
as expected, the mathematics indicates that, in the absence of out- 
flow from the chamber, all of the cardiac output accumulates in the 
chamber. A similar situation exists for V (¢) as K approaches zero. 
However, P(t) approaches P(0) as K approaches zero. As R be- 
comes less and less, on the other hand, 


n 


Via)lpso = napa Ee — na) sin {= - na)} + 


u(t re =) sin {™ (1 na F =)tI} | (17) 


Since y — 0 as R — 0, both the pressure and volume tend to vary 
synchronously with W(¢), In the theoretical situation of no periph- 
eral resistance whatsoever, both V(¢) and P(t) become zero, and 
the cardiac output goes immediately out of the chamber. A similar 
situation exists for V(¢) as K approaches infinity. The pressure 
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has a distinctly finite value as kK becomes infinitely large, the co- 
efficient of the sine terms becoming 6f. The significance of this 
result is not clear. The mathematical results of variation in a, m, 
and 6 correlate well with the expected patterns in the physical 
realm. 

Use of the Laplace transform leads to a representation of aortic 
pressure qualitatively very similar to experimentally recorded pat- 
terns (Figure 3). In addition, it yields equations (12), (14), and 
(15), which may be useful in a practical way. 
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SOME EXTENSIONS OF THE USE OF MATRICES IN 
POPULATION THEORY 


M. H. WILLIAMSON 
DEPARTMENT OF ZOOLOGY, UNIVERSITY OF OXFORD 


(PRESENT ADDRESS: SCOTTISH MARINE BIOLOGICAL ASSOCIATION, 
OCEANOGRAPHIC LABORATORY, EDINBURGH) 


This paper extends Leslie’s vector and matrix treatment of popula- 
tions. A simple matrix is given for species in which adult mortality and 
fertility are independent of age, but in which the juvenile mortality rate 
differs from the adult. The population vector can be changed into apopu- 
lation matrix. This should allow equations using functions of the size of 
the population to be developed. Genetic variables such as sex or other 
polymorphisms can be introduced, and the notation allows different sys- 
tems of selection or non-random mating to be specified. 


Matrices were introduced independently into the study of popula- 
tions by E.G. Lewis (1942) and P. H. Leslie (1945, 1948). Some 
simple extensions of their theory are presented in this paper. These 
are: a smaller matrix than can sometimes be used for premultiplying 
the population vector, the replacement of the population vector by 
a population matrix, and the introducticn of variables other than 
age into the description of the population. Leslie wrote the num- 
bers in the successive age classes of a population of females as 
a column vector n= (no, «++, 7,)° and premultiplied this bya 
square matrix 


Io ii oe iS 
0 
= M, (1) 
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in which the f, are, in a strict sense defined by Leslie, the fertil- 
ity of the various age classes, and the p, are the probabilities that 
an individual in the age class z will live to be in the age class 
+1. Clearly f,20,1>p,;>0. Leslie was particularly interested 
in the way a population approached the stable age distribution 
from an arbitrary age distribution and in the effect of arbitrary dis- 
tributions on the rate of increase. In his second paper he con- 
sidered the problem of limited population growth of the logistic 
type, affecting either mortality or fertility by a factor independent 
of age but dependent on the total population size, and the problem 
of a predator-prey system. 

The matrix M was treated by Leslie as an approximation to the 
true continuous distributions of fertility and death. However, for 
many species living in places where there is a yearly cycle of 
weather or for species which change suddenly from one form to 
another, such as insects, the matrix may well be a better, as well 
as a more tractable representation than any simple continuous 
operator. In some cases the time spent in the young stages is 
much shorter than that in the old.’ Many birds are ‘‘juveniles”’ for 
one year and ‘‘adults’”’ for many, and the death rate and fertility of 
the adults are much the same however old they are (Capildeo and 
Haldane, 1954). Such species can be represented by a system of 


this sort 
Bee n, a f Mo (2) 
Pj Pee Na Pn, + Pan, 


in which the left-hand term is the same as an infinite Leslie ma- 
trix with all p, and f; above p, and f, equal. As the work in com- 
puting a matrix increases rapidly with its size, this variation 
might well be worth while even when the assumptions behind it 
are not really true. 

In elementary theories of population dynamics, that is, those in 
which age or other distributions are not considered, it is natural 
to deal with functions of the population size. Polynomial func- 
tions of vectors do not exist; polynomial functions of square ma- 
trices do, provided that the matrix is not singular, that is, provided 
that its determinant is not zero. It is no use making the column 
vector n into a square matrix simply by repeating the column & + 1 
times, or by using multiples of the column, because those matrices 
are singular. It can be made a matrix N by making it, asit were, 


ie ww ae 
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a diagonal vector, a square matrix with the values of the vector 
along its leading diagonal and zeros everywhere else. This matrix 
will lose its form on premultiplication by M, and after & multipli- 
cations the vector will in many cases be represented by a matrix 
with no zeros. In this matrix the sum of any row gives the total in 
that age class, while each column will give the individuals derived 
from the age class which was in that column in the first ‘‘diagonal 
vector.”’ This will be clearer with an example. Consider a 3 x 3 
M with all f;, p; > 0 acting on a 3 x 3 “‘diagonal vector’’ with all 
n,>0. Then the N matrices in the second and third generations 
are 


for f:m% fas fom. + fzPo%o fofits + f2Pi% folate 
Po 0 0 and PofoMe Poli™ Pola ° 
0 pn, 0 P1Po% 0 0 


Both the matrix M and the vector n refer only to females, but both 
sexes should be included if the mathematical description of the 
population is to be complete. This can he done if the p,, f;, and 
n, are themselves matrices P,, F,, N, with different elements for 
the two sexes. It can be shown that if the 7, are split into row 
vectors (n,3, 7,9), then the p; must be scalar if the Leslie form is 
to be kept; but if the n, become column vectors, then the p; can be 
changed into matrices of the form 


Pt f) 
0 P:9 


PS 


t 


In the simple case of all the F;, except F, being null, it is easy 
to see that the form of the F, matrix will depend on the mating 
system. To take extreme cases, a species in which one male can 
fertilize any number of females will have a different F, than one 
in which the females are monogamous. If there is more than one 
non-null F, their form may well he awkward in calculation as their 
components may be functions of N. (j #i) as well as of N.. The 
F, will anyhow be special cases of the genetic matrices which are 
discussed below. As an example of the use of P,, F,;, and N,, con- 
sider a population of birds on an island with the population system 
described by equation (2). Suppose that the birds are monogamous 
and hold territories of which there are a maximum of T on the 
island, and that the pairs produce equal numbers of males and fe- 
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males. Then the left-hand side of (2) can be enlarged to 

0 0 p o) 53 

0. . Oe 50 ae N59 

P jd 0 Pad 0 Nad 

0 P;2 0 PaQ 1,9 
where ¢ =min (",3, 149 T)f, and where each pair produces 27 
males (regarded as f from the father and f from the mother) and 2f 
females. 

The sex differences are a special case of genetic differences, 
so, as before, the introduction of genetic variables into the column 
vector n will change the n, into column vectors, while the p,; and f; 
become matrices. Evidently the size of both vectors and matrices 
will become impossibly large if more than a few genetic factors 
are brought in, but genetic matrices might be worth using in cases 
of genetic polymorphism where a few genes are responsible for 
most of the visible variation.’ Consider a single gene with but two 
alleles and write d, h, r as the frequencies of the homozygous 
dominants, the heterozygotes, and the homozygous recessives, 
with d+h+r=1. Then the P, will again be diagonal matrices. 
In the simplest case of random mating and all F; other than F, null, 
the genetic part of F, is 


dele Dptdaeke 0 
r+dh gd+dh+ir d+thi. 
0 ar+ch rt+ah 


The elements of the matrix are derived from the expected frequen- 
cies of mating, so this matrix can be adapted to any specified 
form of non-random mating. The rows give the contributions to d, 
h, and r in the next generation from, in the first column, d parents, 
and so on, For instance, the bottom right-hand element gives the 
frequencies with which 7 offspring will be produced from r parents, 
Such parents will, under random mating, mate with other r with fre- 
quency 7 and this will give only r offspring; r will mate with A with 
frequency / and this will give jr and +A offspring. Different fer- 
tilities in the different matings can easily be given in these F. 


If the frequencies of the three genotypes are written d, r and 2A, 
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and remembering that they sum to one, the matrix takes the rather 
simple form 
d+h 3(d+h) 0 
r+h 3 d+h 
0 $(r +h) rth 


> 


or an even simpler one, with p for d+ h and q for r+ h, where p 
and q are the frequencies of the two alleles and p + q = Is 
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A generalization of Landahl’s approximation method (H. D. Landahl, 
Bull. Math. Biophysics, 15, 49-61, 1953) for non-linear diffusion prob- 
lems is suggested. The method is applied to sorption, desorption, and 
free diffusion problems involving concentration-dependent diffusion co- 
efficients. With some limitations, the results compare favorably with 
those obtained by numerical methods. 


Particle interactions often introduce non-linearities into dif- 
fusion systems. Frequently these interactions affect the activity 
and/or mobility of the diffusing species and as such the non- 
linearity becomes manifest as a concentration-dependent diffusion 
coefficient. More complicated non-linearities can occur in systems 
where there are several interacting species as, for example, in an 
electrolyte mixture. Analytical solutions to these problems are 
known only in a few relatively simple cases. 

An approximation method for transient diffusion where the dif- 
fusion coefficient D remains constant has been recently proposed 
by Landahl (1953). Briefly, this method can be illustrated by con- 
sidering the sorption problem in a semi-infinite medium, that is, 
with the boundary conditions: 


c(z,0) =9, 22 0; (1) 
c(0, t) = os t? 0; 


where c, « and ¢ denote the concentration, space coordinate, and 
time, respectively. The method then consists of assuming a linear 
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concentration profile (see Figure 1) and determining the position 
of a time-dependent x intercept r(¢) by applying Fick’s law to- 
gether with the mass continuity principle at the boundary 2 =0. 
Landahl applied this method to the above and to other boundary 
value problems and obtained results which compared favorably 
with corresponding exact solutions when they were available. 


Co 


X=0 X=r (t) 
Xx 


FIGURE 1. Sorption. 


If this method is applied directly to a case involving a concen- 
tration-dependent diffusion coefficient, it will be found that the 
solutions do not indicate this dependence. Further, the fact that a 
linear concentration profile is a good approximation in the case 
governed by a linear differential equation in no way assures us 
that it will be adequate in the case described by a non-linear 
equation. However, upon re-examining the method, it can be seen 
that if D is constant, a linear concentration profile is simply the 
steady-state concentration distribution between z= 0 and x = r(z); 
and hence the method could be thought of essentially as a quasi- 
Steady-state approximation. In other words, from this point of 
view it is assumed that as equilibrium is approached, the system 
passes through a sequence of steady states. The particular steady 
State that exists at any moment is determined by the mass con- 
tinuity condition at the boundaries. Reformulated in this way the 
method may be applicable to both linear and non-linear problems. 

In this paper we shall compare solutions obtained by the quasi- 
Steady-state method mentioned above with some available numeri- 
cal solutions. We shall restrict our attention to problems arising 
from the simple Fick equation given by the expression: 


a 
Je tt tayee (2) 
Ox 
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where J is the net flux across a unit surface perpendicular to the 
a axis, and D(c) is the diffusion coefficient. We reserve for a 
subsequent paper the problem of diffusion in a mixture of more 
than two ions and consider here only those non-linearities that 
can be entirely incorporated into a concentration dependent D. 

The method will be applied to three related boundary value 
problems. [Many numerical solutions to these can be found in 
Crank (1956)]. For other types of boundary conditions variations 
of the method analogous to those used by Landahl for correspond- 
ing cases probably could be employed. First we shall derive the 
results for these problems and then compare them with numerical 
solutions. 


Sorption. Consider the sorption problem illustrated in Figure 1. 
Here the boundary conditions are given by (1), but in addition, the 
approximation also yields the relation: 


c(a, t) = 9, 22 r(t). (3) 


If at time ¢ a steady state exists between x = 0 and 2 = r(¢), then 
J is independent of « within that interval, and (2) can be written 
in the integrated form 


Ja= [-"D(Ad0= Dol (c), (4) 


where the function /(c) is defined by the identity in the above ex- 
pression, and D, = D(0) is the limiting diffusion coefficient at 
infinite dilution. 

When « = r(é), then c = 0 and we have from (4) 


ee 1(0) Dy 


i 


(5) 


The continuity condition at the surface x = 0 can be expressed by 
equating at each moment the cumulative amount of material which 
has diffused across this surface to the amount of material that 
exists beyond it. That is, 


. r 0d 
[see f ode = | fae de (6) 
dc 
0 0 co 
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Using (2) to calculate da/dc and (5) to replace J, and then dif- 
ferentiating (6) we find 


wae = row [oer ao < 10) (7) 
Solving (7) for r(¢) with the condition that 7(0) = 0 we have 
r=2VKt, (8) 
where 
spt D,1(0)? (9) 


2 Co D(od 
—= c Cc) ac 
Dy 0 


If D(c) = D’ where D’ is a constant, then K should and does reduce 
to D’, the same solution as that given by Landahl. 

Further discussion will be expedited if we introduce the variable 
y by the definition 


oo 


2/D,t’ 


y= (10) 


and the function M by the definition 


U(c) = ra [encoae] oo (11) 
00 


In these terms the concentration c at any 2 and ¢ can be written as 
a function of the single variable y. Combining (4), (5), (8), and 
(9), c can be written as the root of 


Ic) 5 er ee 
4 (c,) y ’ ¥~ OVD,’ 
12 
te (12) 


c=0 2 . 
scgecaa Te 


Desorption. If we consider the corresponding problem of desorp- 
tion in a semi-infinite medium (see Figure 2), then, instead of (1) 
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and (3), we now have 


c(x,0) = ¢5, 
c(0, 2) ve ; (18) 
and 
e(a, t) = cy, 22 r(t). (14) 


Proceeding, with these modifications, in the same manner as 
above, we arrive at the following expression relating c to y: 


pe OHH yg 
Dc, 1(0) —M (co)? SS Nee 
(15) 
= r(t) 
= Co? y= 2/D,t 


Free diffusion in an infinite medium. Finally, we may consider 
the problem specified by the boundary conditions 


e(z, 0) = co; x < 0; 


c(z,0) = 0, z>0. (16) 


Referring to Figure 3, we now define the points x=r(t) and 
2 = —s(t) by the relations 


c(z,t)=c ’ 2 <— s(t); 
c(z, t) = 0, a> #(t). } (77) 


Note that c, no longer represents a fixed concentration at the 
ordinate 2=0. In addition, the approximation stipulates that, 
within the interval —s(¢) <z<-r(t), J is independent of a. Then, 
within this interval, we may again define the function /(c) by the 


ot 


X=0 Xe=r (t) 
Xx 


FIGURE 2. Desorption: FIGURE 3. Free Diffusion. 
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relation 


Has s)= [-°D(0)dd = Dye) (18) 


Let c,’= ¢(0, 2), i.¢.,; cis the concentration at 2 =0.: Then the 
expressions 
Fee T(9) — gs é Meo) p 
r s 


(19) 


0? 
from which 
r I(0)- ites) 


< eto ae (20) 


are obtained from (17) and (18) by letting e =r andz =0. 
The continuity condition at z = 0 requires that 


t r 0 
[ sae= [ ede = [ (c, — c) dx. (21) 
0 0 -s 


Using these relations and proceeding in the same manner as in the 
previous cases the expressions 


. edie = 
Meg)? BD YB YDS? 
a, = r(t) : 
c 12 a (22) 
- (2) 
C= Cp) ys aD = 
0 
with 
2 pes 7 
M(c,’) = ee oP (c) do ’ (23) 
Oo 


are readily obtained. 
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It remains to determine the value of c,. For this purpose, the 
last equality in (21) can be rearranged in the form 


i = la S ie ar) 0 
¢~= Cc —_dac?:er—oOrororr————————— = 
: . | de DylN(0) Mey) | ae hay 


s 


0 
cy | dz = Cy & (24) 


Ta} 


Comparing (24) with (20), we find that c,’ can be determined from 
the relation 


2 
ee ie. (25) 


and thus our free diffusion approximation is completely determined 
by (22) and (25). 


Some examples and their comparison with numerical solutions. 
If the dependence of D on c is known explicitly, then the function 
I(c) is readily computed from (4) and M can be obtained from (11); 
Knowing / and M, the solution to the sorption problem is given by 
(12), the solution to the desorption problem by (15), and the solu- 
tion to the free diffusion problem by (22) and (25). 

In Table 1 we list four different forms for D(c) and the resultant 
expressions for / and M, all in terms of cy and an arbitrary param- 
eter qa. These are then used to calculate the dependence of c on y, 
and the results are plotted in the figures (solid curves), where 
they are compared with numerical results compiled by Crank (loc. 
cit.; broken curves). 

In making this comparison, it should be remembered that the 
approximation assumes that the concentration distribution can be 
truncated. Thus, for example in the case of sorption, instead of 
the stringent but approximate expression c(z, ¢) = 0 for x = r(¢), the 
exact solution requires only that 

lim c(a, t) = 0. 

x90 
Consequently it should not be surprising to find that the method 
becomes more and more inapplicable as z — r(¢), and indeed this 
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is revealed in the figures. However, the figures also reveal that 
outside of this region the numerical and approximate solutions 
closely agree, and that similar remarks could be made for desorp- 
tion and free diffusion. 

All of the examples illustrated in Figures 4A, 4B, 5A, 5B, and 
6 are for positive values of the parameter a, and this means that in 
all these cases Dis a monotonic increasing function of c. In 
order to establish the method, it would be desirable to consider 
decreasing functions of c. However, as pointed out by Crank, this 
is always implicitly accomplished if both sorption and desorption 
have been considered. This follows from the fact that if (c) — c) 
is written for c, the diffusion equation remains unchanged but con- 
ditions (1) change to conditions (13). Consequently the exact 
sorption solution when D is a given function of c is also the exact 
desorption solution when D is that same function of (c, — ¢) and 
vice versa. It can be readily verified that the same substitution 
transforms the approximate sorption solution into the correspond- 
ing desorption solution and vice versa. These remarks justify the 
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FIGURE 4A. Sorption. D= D) ( +a 3 . The numerical solution is 
a 
“0 
represented by the broken line. The approximate solution is represented 
by the solid line. 
c 
FicureE 4B. Desorption. D=D) {1+a@—}- The numerical solution 
C 
0) 
is represented by the broken line. The approximate solution is repre- 
sented by the solid line. 
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FIGURE 5A. Sorption. D =D, exp (2 a . The numerical solution is 
represented by the broken line. The SG solution is represented 
by the solid line. 

FIGURE 5B. Desorption. D =D, exp f a . The numerical solution 
is represented by the broken line. The ee oe solution is repre- 


sented by the solid line. 
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sented by the broken line. The approximate solution is represented by 
the solid line. 
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extension of our general results (for sorption and desorption) to 
problems involving diffusion coefficients which decrease with 
concentration. 

In Figure 7 a single example of the sorption problem is illus- 
trated where D first increases to a maximum value approximately 
seventeen times its initial value and then decreases back almost 
to its initial value again. 

Finally, in Figure 8 the problem of free diffusion is illustrated 
for coefficients which increase with c and for one which decreases 
with c. 

Some idea of the generality of the method is obtained by an 
examination of the extent of departure of D from constancy in these 
examples. This is indicated for each curve in the figures by the 
value of k, here defined as the ratio of the greatest value of D to 
its least value within the concentration range under consideration. 
Thus, for example, in Figure 4A with a = 25, D ranges from D) to 
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cal solution is represented by the broken line. The approximate solution 


is represented by the solid line. 
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FicurE 8. Free Diffusion. D= Dp) exp | @ =| . The numerical solu- 


a) 


tion is represented by the broken line. The ‘approximate solution is 
represented by the solid line. 


26 D,, and accordingly & = 26. It will be seen that values of & as 
high as 200 have been included in the examples. 

Despite the diversity of types of diffusion coefficients studied 
and despite the strong departure of D from constancy in many of 
the examples, the results are essentially similar; and hence the 
accuracy of the method appears to be reasonably independent of 
the particular form taken on by the function D(c). Considering the 
simplicity of the method, the results are perhaps surprisingly 
accurate. It is, of course, limited in that it cannot predict the 
dependence of c on z in the neighborhood of the truncations 
[ew = r(t), « =—s(t)]. However, for many problems, these values of 
c play a rather insignificant role. For example, in the sorption 
problem consider the cumulative amount of material that has dif- 
fused through the surface z=0 up to time ¢ If we denote this 
quantity by Q, we have 


Q -/ cdz = 2/D,t i cdy. (26) 
0 0 


Apparently then, the exact method and the quasi-steady-state ap- 
proximation differ only in the value they assign to the integral on 
the right-hand side of (17). For the quasi-steady-state method, the 
value of this integral can be computed directly and is simply equal 
to (1/2)M. Thus our approximation method yields Q =MV/D,¢é. 


For the cases considered, values of aa obtained by the ap- 


proximation are compared in Table 2 with those obtained with the 


Bo 
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TABLE 2 
Q Q 
2D, t 2V Dt 
D(e)/Dy a APproxe exacte % Error 
2.5 0.82 0.92 10.9 
rae Ged 1.23 1.36 9.6 
1+ a— 25.0 2.05 2. 30 10.9 
i) 50.0 2.92 3.20 8.8 
100.0 4,12 4.51 8.5 
0.69 0.63 0.71 11.3 
: 1.61 0.89 0.98 9.2 
a. Co 2.30 1.15 120 8.0 
3.91 2.18 2.33 6.4 
4.55 2.86 3.03 5.6 
5.30 3.93 4,13 4.8 
: 0.80 0.80 0.85 5.9 
0.90 0.93 1.02 8.8 
- 0.95 1.07 1.13 5.3 
1- a— 0.98 1.23 1.33 (es) 
0 
1 
1.32 1.43 8.4 


a) % 


Cc ye, 2 
1 — 3.29 — + 2.88 si 
help of the numerical solutions compiled by Crank. It can be seen 
that the quasi-steady-state approximation consistently underesti- 
mates Q (as might be expected by the error introduced in truncating 
the distribution). Further, the average error is only about 8% and at 
most is slightly greater than 11%. If this magnitude of error is not 
tolerable it can be reduced by noting that the method always under- 
estimates and that the range of error is between 5 and 11 per cent. 
Hence, if an empirical correction of 8% of the predicted value of @ 
is added on to Q, the error probably would be reduced to within 
+4%. 

If the problem requires more accurate results, it is suggested 
that the quasi-steady-state solution discussed in this paper would 
form a desirable first approximation to be used in iteration proce- 
dures. (cf., Crank, loc. cit.) 


I would like to thank Dr. Clifford S. Patlak for reading and dis- 
cussing this paper. 
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Recently a theorem for representing current generators in a volume con- 
ductor by the superposition of a central dipole, quadrupole, octopole, 
etc., has been established by G. C. K. Yeh, J. Martinek, and H. de Beau- 
mont (Bull. Math. Biophysics, 20, 203-16, 1958). This theorem makes 
possible the representation of any discrete or line, surface- or volume- 
distributed current source by a unique model which can be determined for 
each given case by surface potential measurements and closed form 
analysis. In this paper the multipole representations of an eccentric di- 
pole and an eccentric double-layer are obtained in terms of the various 
parameters of the assumed singularities, and the contributions to surface 
potentials due to each of the multipoles are compared. Certain numerical 
results corresponding to those of E. Frank (Amer. Heart J., 46, 364~78, 
1953) are carried out and compared. Furthermore, the multipole representa- 
tion of a partially damaged double-layer is also determined and compared 
with that of an undamaged one. It is concluded that within the range of 
parameters corresponding to human subjects the higher-order multipoles 
can contribute significantly to the surface potentials compared with the 


dipole. 


Multipole expansion of any potential in an infinite domain. -In 
Yeh, Martinek, and de Beaumont (1958), it has been shown that the 
potential due to an arbitrary distribution of current source density 
Bio, Yor Zo) which is in the domain 7 < a, at a point outside this 
domain can be expressed by the following series: 


0 (cy) (z) 
$y = 1 98> + D, 5° + Dy bo" + 92 bo + 
( (zz) 
Ox oot + Oy ie z Q, z Po Zs 
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In equation (1) ¢ is defined as 


q ={[[ode,dyodeo, 


(1) 


(2) 


and is the strength of a source at the origin whose unit potential is 


1 


Tr 


(0) _ 
0 


D,, D,, and D, are defined as 


dD. = | [feodes dy, dey, 


v 


d= |[[yspdeydy, dey, 


v 


D, = [[[eoedey dy, dey, 


v 


(3) 


and are the strengths of the three components of a central dipole 


whose unit component potentials are 


EO) 

0 On * ae? 
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Qin> Vyyr+++Q,,.++- are defined as 


Qex = |[[=2 eae, dy, dz, 


i (5) 
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and are the strengths of the nine components of a central symmetric 
quadrupole whose unit component potentials are 
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and are the strengths of the 27 components of a central symmetric 
octopole whose unit component potentials are 
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Ge 0-6 6 bie 8 6 elie: ©: (@ feo 0, 0 ie 07.0 Je 6; 68 0) 0 6: 10 'o 16 e810) 6 le eee a ean eae eee ere 


ig 


Shi Ne 
ee = ayes as ep = ie Ciel = Patek: z Pace = 
1 g° (=) 5axy2 


nt, dxdydz rt] Ort s 


The succeeding terms in equation (1) correspond to potentials due 
to a central symmetric 2*-pole, a central symmetric 2°-pole, etc. 


Multtpole representation of an eccentric dipole. If the current 
source density p(%), ¥), 2,) is such that a single sink of strength 
~ 2 is located at the point (0, 0, 5) and a single source of strength 

é 
+ Ee is located at the point (0, 0, 5+.) and if we let e« become 

€ 
vanishingly small, we have an eccentric radial dipole of strength 
#., pointing towards the positive 2 axis at the point (0, 0, 4). To 
represent this dipole by a combination of central multipoles we pro- 
ceed to determine the corresponding strengths of the components of 
the multipoles here. From equation (2) the strength of the central 
source is 


q=lim (ts + Es) ~0, (9) 


€& 0 e € 


From equations (3) the strengths of the dipole components are 


rn r 
D,, = lim - tta,+ He,| =0, 
Xx,=0 


€>0 € € 
; lL le 
D. =) <i Uz og 
yotin [- 90+ ae 0, (10) 
= By Ht 
3, = lim [20+ Stee a] mae 
0 


Sr a ee e ANA 


NN es 
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From equations (5) the strengths of the quadrupole components are 


ee a awe 2,2 = 0» (11) 
Q,, = lim | Sad + (2, + 07| =2bu 
zZou=mb 


yYvY xxY xyy xz “ZZ yzz yz “YZ 
: LU 
O inh er a a (8 ee)” = 3)? 
Zzz €->0 0 € ( 0 ) eis Pz» 
0 


and so on, for higher-order multipoles. 

Thus the dipole strength D, =, indicates the total moment of 
the eccentric dipole, the quadrupole strength Q,, = 2p, indicates 
the eccentricity of the eccentric dipole, and the octopole strength 
O22 =? b? 1, indicates the square of the eccentricity, etc. 
According to equations (1), (4), (6), and (8), the potential due to 
this eccentric dipole in an infinite domain can be expressed as the 
sum of the potentials due to the central multipoles with component 


strengths obtained in equations (9) to (12) as follows: 


cos 6 3 cos? 0-1 
Boat a OS ley k a 
cos 6(5 cos” 6-3) Ly b 
Ua er era pe ee cos 6+ 3 (8cos? 6-1) (18) 


8° LB : 
roo cos (5 cos? 6-3) +... for 74> bs 
. 


a 
where the relations 
pimrisin Uicos @, 94 =t-ain O'sin Y 5. 2= 7 COS 0, 


have been introduced (r, 6, ~ are spherical coordinates). It is also 
known (see, for example, Yeh, Martinek, and Ludford, 1955) that 
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the exact potential due to such an eccentric dipole is 


(cos 6 - | 
uw, (2 — 4) H J 


dp = ao a” te eee 
a3 [w? + y? + (2 - 5)? ]? : 1-2 ~) eos o-(-F 
r 


b 
Thus at the point 6=0, r= rere equations (13) and (14) show that 


the central dipole, central quadrupole, and central octopole con- 
tribute 36%, 28.8%, and 17.28% respectively of the total potential 
due to the eccentric dipole along 6=0 at r= 5, whereas at the 


b 
point 6=0, r= 09 they contribute 64%, 25.6% and 7.68% respec - 


tively of the total potential. At the points r — ~ the central dipole 
alone can account for the total potential due to the eccentric dipole. 

Let us now consider an eccentric dipole of strength p, parallel 
to and pointing towards the positive z-axis at the point (0, 0, 3). 
This is not a radial dipole and can be obtained by locating a single 


x 


sink of strength — RP the point (0, 0, 5) and a single source of 
€ 


x 


strength + r= at the point (e, 0, 5) and letting « become vanishingly 
€ 


small, The component strengths of central multipoles representing 
this dipole can be determined by a method similar to that used for 
the radial dipole demonstrated above. From equation (2), 


q = lim (- &s + ts) m0. 
€>0 € € 


From equations (8), 


‘ u a 
D. = ij Fae Li = 
é lin [ = Ont S(t of Bos 
0 
: a UL 
y eat : € Yo - € vo) =O a 4h) 
0 


€>0 


D, = tin |- #5 24+ #2.) =0, 
€ 
Zo™b 
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From equations (5), 


eee 
qo2 8 |- i+ te, +9) =O, 
X%9=0 
Qyy = 922 = 2ey = Qy2 = 9; (16) 
oC ae eed Px 7, 
@,, > lin | : 22> + : (a +926] = by. 
%o™0 
Zo™u™b 
From equations (7), 
taf = 2 |- = = + “£(@ + a’ =0, 
%o™0 
Seve Jae = Le = Lae im Ores vez LP = Ons =0, (17) 
m lL 
Ole = lim |- aed Hs (ey + 029 | = b7n,, 
X%o™= 0 
Zo™m5 


and so on, for higher-order multipoles. 
Hence from equations (1), (4), (6), and (8) the potential due to 
this dipole in an infinite domain can be expressed as 


sin 6 cos 5 3 cos § sin 6 cos 
——————— - ss 


Pop at ie 9 a Hy 3 


x rT a 

> 3 sin 6 cos (5 cos” 6 - 1) 

ee 
= ort 


: 5 (18) 
—~ sin gcos ¢|1+8(~)cos 0+ 
: r 
a/b) 
sae (5 cos? 6-1)+...| for r>d, 
* 


where the first, second, and third terms are potentials due to the 
central dipole, central quadrupole, and central octopole respectively. 
It is also known (Yeh e¢ al., 1955) that the exact potential due to 


such an eccentric dipole is 
HB, % 


ote en as 
pes bee ey tt (a - 8)°L 


Po 
(19) 


ie 1 
— sin 6 cos 
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Thus along the line 6 = 0 there is no potential due to this dipole. 
b b 

. pe if bee Ps eS 

On the circle @= cos a r A 
quadrupole, and central octopole contribute 77%, 36.96%, and 
-3.7% respectively of the total potential due to this eccentric di- 


b b ; 
pole. On the circle 6 = cos ! Pe r= ao they contribute 94.06%, 


11.29%, and -4.51% respectively of the total potential. And at the 
points 7 —> « the central dipole alone can again account for the 
total potential due to the eccentric dipole. 

Similarly an eccentric dipole of strength y, parallel to and point- 
ing towards the positive y-axis at the point (0, 0, 6) can be repre- 
sented by the combination of a central dipole of component 
strengths 


, the central dipole, central 


Di=p,, D, =D, =9; (20) 
a central quadrupole of component strengths 


Ce ee b py Qe Pai Ose = Ae = Sie ? Q.2 ca 0; (21) 
a central octopole of -omponent strengths 


Ove, = Cmte O 


KX OTe Py La = eo =O = 0 rs 


Ve x¥¥ xx Zz (22) 
0 zz a = Oe =O; 


x 


and so on, for higher-order multipoles. The potential due to this 
dipole in an infinite domain can be expressed as 


p ee cinta aid peal Veaee 
OP rag? ~ ' 


Pere (23) 
> (2) (5 cos? é-1)+.:.| for r> 6, 
2i\F 
as compared with the exact potential 
boom ; 1 
Pou, = =a Sin O sin yp: (24) 


Lo. Ae aes 
- 2|—]cos 6+|— 
r r 
Finally, a general eccentric dipole at the point (0, 0, 4) with 
component strengths y.,, ty, and w, can also be represented by a 


combination of central multipoles whose component strengths are 


given by superposition of those for the three eccentric dipoles dis- 
cussed above. 


Ne Ree ee Ae a OAT = Gee 
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Multipole representation of an eccentric double-layer. Let us 
consider a double-layer current source of uniform moment M per 
unit area on the spherical surface r = A with a circular rim of radius 
c whose center is at (0, 0, 5) and whose axis coincides with the 
z-axis as shown in Figure 7 of Frank (1958). We may obtain this 


M 
double-layer by placing a simple layer of uniform strength — — per 
€ 


unit area on the spherical surface r = A with a circular rim of radius 


e in the plane z = 6 and by placing a simple layer of uniform 
strength 


M’ M h(h — b) 


gh > g ° Sy SA Pees ear A 
(A+ e)(h+e-b-€- 
h 
per unit area on the spherical surface r = A+ with a circular rim 


of radius c ( ~ +) in the plane z = 6 ( - a and by letting « be- 


come vanishingly small. 
Using the spherical coordinates we have from equation (2) 


ampecos tb/h M 
ffm Lim [ | — —h? sin 6d0dq + 
€ 
0 


€>0 
0 
2m7pcos + b/h M 
Fe if if Sees (25) 
€>0 0 Yo € 


hoa) ne (h + 6)? sin 6d6dy =0. 


(he (n+e-b-e5] 


From equations (3) we have 


2mrcos 1b/h M : 
D, = tim | | h sin acoso (- *)i sin 0d 6d + 
i€ 
0 vo 


€>0 


2m7pcos 1b/h M’ 
lim [ E: (A+ €) sin 6 cos (n+ 9? sin 6d 0d = 
€ 
0 ~0 


€>0 


(26) 
1 3 
lim — [- MA° + 


€>0 € 


cos !b/h 27 
M’(h+ olf sin? ao | cos pdg=0, 


0 0 
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D,=0, 

2mpcos 1b/h —M 
p,= tim | [ h cos 0|\— JA? sin 0dOd@ + 
Ea eA el € 


2mpcos 1b/h M’ : 
lim if [ (i+ €) cos 6 @G +6)? sin 6d6dQ = 
€>0 0 € 


0 (26) 
(cont.) 
bth eh 
2M lim [ Sl atth: 
€>0 1 € b 
(1 +e-b-e z 
h 
cos? 6 eo 
cos 6d cos 0 = — 27h*M = 7M (h? - 6?) =a7Mc?. 
1 
From equations (5) we have 
27 pcos 1b/h —M 
Qe Lim [ il h? sin? @ cos? 9| — }A? sin 6dOdQ + 
€>0 0 0 € 
: 27 pcos 1b/h M’ 
lin [ [ (1+)? sin? 6 cos? p(—](A+ 6)? sin 6d6dQ = 
€>0 0 0 € 
-h* h(h-b 
M lim Steep eee nee is) 
e690 | € b 
€ (i +e-b-e 4 
h 
(27) 


an cos 1b/h 
i cos? eae f sin? @ sin 6d6 = 
0 0 


2h3 ~ b(3h2 — b?) 
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b/h 
M-2A8 nf (cos? 6 —1)d cos 6 = 2nMh8 
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2 
5 7M [2h ~ (3A? — 8%), 


2 
Qyy= y TMIZA® - BBA? - bY], 
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: 27 pcos 1b/h —M 
= Lim [ ip h2 cos > o(— ye sin 6d6dq + 
0 
7 pcos }b/h M’ 
Sap if (A + €)? cos? (= ~)¢ +)’ sin dd0d9 = 
0 


b/h 
2m -2A%. f cos? 6(—d cos 6) = 
1 


Apt 4 
47Mh8 . Sone. 
7 a 5 7M (a* - 4), 


: 27 AO Eoe : M ae 
Q. = Lim / h* sin* 6 cos @ sin (=) sin 6d 6d + (27) 


(cont.) 
27prcos 1b/h M’ 
Lim [ ib (A + €)? sin? 6 cos 9 sin o( )x 
€>0 é 


0 0 


(h +)? sin Odédy =0, 


277 pcos *b/h 
Q. -= lim [ ih h? sin 6 cos 6 cos 9 | — ree sin 0d Ody + 


€>0 


27 pcos 1b/h M’ 
lim [ | (h + 6)? sin 6 cos 6 cos 9 =) x 
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From equations (7) we have 


27 pcos 1b/h : —M E 
Oven = Lim [ [ h? sin® 6 cos® 9|— ]A? sin 6dddy + 
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€>0 
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0 Yo € 


=k Ath SBA : 
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Qm7pcos tb/h —M ; 
Omi him [ h® cos? 6| — JA? sindddd® + 
Z22Z €>06 “ 4 € 


2 cos !b/h M’ 
lim [ ii (h+«)* cos? (aso sin 0d6dQ = 
(3 
o 0 
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€ 
0 
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E90), Jo € 
27 pcos” 1b/h M’ 
lim { [ (A + )° sin? 6 cos 6 cos? (=) x 
€>0 0 0 € 


(h +e)? sin 0d6dQ = 
277 b/h 
M-3h4 -f cos? yd ‘1 sin? @ cos 6(—d cos @) = 
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= 7 rm (a? ~ b)9, 
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= 77M (A? ~ d9)9, 


ECCENTRIC DIPOLE AND DOUBLE-LAYER 45 


Ony27 9: 
: 27 pcos !b/h —M 
0. ..= lim { i ® sin 6 cos? 6 cosy ie x 
ero), Jo € 


h? sin Od6d@ + 


h 2amrcos lb/h M? 
Kim | il (A +)? sin 6 cos? 6 cos o(“] x (28) 
0 Yo 
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(h + €)? sin Oddy = 


277 cos 1b/h 
Mo3nt | cos pa | sin? 6. cos? 6d6=0, 
0 0 


Thus, as far as the potential outside the domain r=h is con- 
cerned, the eccentric double-layer can be represented by the com- 


bination of a central dipole with a non-zero component D, = 7Mc’, 


: 4 
a central quadrupole with non-zero components Q,, = 5 7M (A? — 53) 
2 
and 9, =4,,= = 7M[2h® — 6(3A? — 6?)], a central octopole with 


non-zero components O,,, = z aM(h*—b*) and O,,,=0,,,= 


3 
aM (h? =67)7., ., ote. 


In a similar manner, we may represent an eccentric double-layer 
disc (of uniform moment M per unit area) which coincides with the 
rim of the given double-layer by the combination of a central dipole 
with a non-zero component D, = 7Mc?, a central quadrupole with a 
non-zero component Q,, = 2 7Ma(h? — 6”), a central octopole with 
non-zero component O,,, =37Mb?(h? - 57) and O,,,=0,,, = 


1 
7 Mn (h? — 07)?..., etc. 


With consideration of the expressions of 4%”, 6°”, 62”, of, 
$6", o7... given in equations (6) and (8) in relation with the 
total potential expression given in equation (1), it can be shown 
that both the double-layer and the double-layer disc produce the 
same potential distributions outside the domain r = h as that due to 
the combination of a central dipole of a non-zero component D, = 


46 G. C. K. YEH AND J. MARTINEK 


7Mc2, a central quadrupole of a non-zero component Q.,=2 m™Mc?b, 
3 
a central octopole of a non-zero component O,,, = gmc" (4b? —c?), 


etc. Thus the dipole strength indicates the total moment of the 
double-layer, the quadrupole strength indicates the eccentricity of 
the double-layer, and the octopole strength indicates something 


about the size of the double-layer, etc. 
Hence the potential due to either the double-layer or the double- 


layer disc can be expressed as 


F cos @ 8 cos? @-1 
doy = Mera 72 sere eae eT ae 


cos 6(5 cos” 6 - 3) | 
—_—_—_————— +... |= 
2r4 


3 
4° be = c*) 


Mre? b 
ne [cos 6 + (“os cos? 6-1) + 
r 


2 {4(2) -(2) } cos 6(5 cos? 6 — 8) + za for r>h. (29) 
8 - r 


It is well known (see, for example, Lamb, 1932, p. 212) that a 
vortex ring is equivalent to a uniform distribution of double sources 
over any surface bounded by it. The axes of the double sources 
must be supposed to be everywhere normal to the surface, and the 
density of the distribution to be equal to the strength of the vortex. 
Hence both the double-layer and the double-layer disc can be re- 
placed by a circular vortex ring of strength M on the rim of the layer. 
The exact potential due to such a vortex ring can be expressed as 
follows (see, for example, Yeh e¢. al., 1955): 


rcos 6-5 
doy = ~ 2M [ee oa (a) + NG X| : (30) 
1 
where 
A(e,k’) = K(k) E(e, k’) + [E(k) — K (Re (31) 
and 
«=F (y, k’). (32) 


In equations (30), (31), and (32) K(k) and E(k) are complete ellip- 
tic integrals of the first kind and the second kind respectively for 
modulus &: EF (e, k’) is the incomplete elliptic integral of the second 
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kind of the argument e and for the modulus k*: F(y, k’) is the in- 
complete elliptic integral of the first kind of the argument y and for 


the modulus k’ and 
ri? = 172 + 2r(c sin 6-5 cos 6) + (62>4 ¢7); 


r,2= 1? —2r(c sin 0 + b cos 6) + (b? +c), 


1 (33) 
pe lre tye Saat 
ei ra 
cos y:sin y:l1=rsin @—c:6b-rcos @:r,, (-7S yn). 
Thus on the line 6 = 0 we have 
2 eg aks 
pat. =(r- 6)? +c7% =r, k=0, k=1, 
y = Er + sin”) — for rz b; 
"9 
eG eae 
e=F(y, 1) = #2K(1) + F \sin : , 1) ~~, 
0 
A (e, 1) = K(O)E (61) ¥2 lim (H (B) ~ K (RL K(R) = 
7 : 7 r—b 
5 F 2E (A) + sin (y tml+0= 5 =e aeeP (34) 
au r—b[a\_ 9) 
Pou = ~ pees er 
r—b 
2M | =| = 
"9 
- b 
Mnc? 2 sa 
7C = 256 
ete es 
Gy He-2) +0) 
r r r 
; b 
For the case of 6 = c and at the point 9@=0, r= TW isis have 
Mac? 2 | 0.6 -1| - 
OM ~ ~ 2 0.4)? 2 234 
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Hence equations (35) and (29) show that the central dipole, the 
central quadrupole, and the central octopole contribute 47.64%, 


38.11%, and 17.15% respectively of the total potential. 
b 


Now for the case of 6=c and at the point 6=0, r= Oe. we 


have from equations (34) 


“Mao? 2 0.8 Mie? 
ou = a on OEE - | = 1.4928 — -. 
r (0.2)° | ./(0.8)2 + (0.2)? 


Hence equations (36) and (29) show that the central dipole, the 
central quadrupole, and the central octopole contribute 66.99%, 
26.80%, and 6.03% respectively of the total potential. 


(36) 


G 
For the case of 6 =0 and at the point 6=0, r= ne we have 


from equations (34) 


= 2 M 2 
(fe ie | ~o.071 te 
r (0.2)° | 12 + (0.2)2 r 


Hence equations (37) and (29) show that the central dipole, the 
central quadrupole and the central octopole contribute 102.99%, 0%, 
and ~ 3.09% respectively of the total potential. 

And equation (29) shows that for points at r —> o the central 
dipole alone can account for the total potential. 

It is interesting to compare the strengths of the multipole repre- 
sentations of this double-layer D, = aMc?, Q..™ 27Mc? b, 0.227 


aro? (4 b? —c*)...etc., with those of the eccentric radial di- 


pole of strength »,, D, =,, Q,, =2bu,, O,,,=367n,...0etc. As 
long as they have the same total strength p, = 7Mc?, they will 
have the same central dipole. If they also have the same eccen- 
tricity 6 they also have the same central quadrupole. But the 
octopole strength for the double-layer is always smaller than that 
for the eccentric dipole unless c — 0, in which case the double- 
layer approaches the eccentric dipole. This shows the importance 
of the higher-order singularities for revealing details about the cur- 
rent generator. To consider the central dipole alone, as has been 
done in the standard leads and vectorcardiography, is consequently 
too rough an estimate. 


ear NY 
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Multipole representation of a partially damaged double-layer. 
Suppose an axi-symmetric portion at the apex of the double-layer 
considered in the last section is damaged and ceases to function. 
The current source then becomes a double-layer of uniform moment 
M per unit area on the portion of the spherical surface r= h be- 
tween the plane 2 = } and another plane z = 6’ > 6. We may obtain 
this partially damaged eccentric double-layer by superimposing a 
double-layer current source of uniform moment ~M per unit area on 
the spherical surface r = A with a circular rim of radius c’ = 
Vh2 — (8’)2 < c whose center is at (0, 0, 6’) and whose axis coin- 
cides with the z-axis over the original double-layer considered. 

The multipole representation of such a damaged double-layer can 
be immediately obtained from equations (26), (27), and (28) as the 
combination of the following non-zero components: 


D, = 7M[c? ~(c’)7], 


2 
Q.. = Qyy = = TM (H ~ BBA? ~ (BY? ~ 5d’ — 871, 


Qe, =e wM{(o)? - B91, 
: (38) 
O = 5 MAES — O41, 


3, tay Ape SOG SSRI NE Se CMON te NE RCS da i ae ot ad 


a ee BS I OEP ETP ONE ONE OT ee Pe DE Naa es Aa 


This system produces the same potential distribution outside the 

domain h =r as that due to the combination of a central dipole of a 

non-zero component D, = 7M (c? - c’?), a central quadrupole of a 
b°2 + bb’ — c? 


Dod 


a central 


12 
non-zero component Q,, = 27M (c? - c’*) ; 


3 ’ vd 
octopole of a non-zero component Oe a pat (c? — c’7)(5 b°? - 


b2 ~ 6c”), etc. Thus suppose c is given. The simultaneous so- 


Jerret- 07 


O 4 
lution of fs and — render } and 6’. Then c’ = 


Zz 
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and the unit strength of the double-layer is immediately given by 
D 


z a(c? — 0°?) 
The potential due to this partially damaged double-layer can be 
expressed as 


ee Ose b°2 4 bb’ ~c? 3 cos? 6-1 
Porm = Male SC, ) 72 a ee Digi eee 
3 cos 0(5 cos” 6 — 3) | 
= (5 5°? — 5? 6 69) «eee 
+ (50 0?) = " 


b’ b \ [b’ c\? (39) 
——— + = — ao — 
r r}\ Pr r 

os 6 +4———[—|_——_— 

3 b’ \2 b \? c\? 

—45|—]) -| —] -6|—] (cos 6(5 cos? 6-3)+... forr>A. 

8 r r ‘a 


The exact potential can be expressed as that due to two vortex 
rings in the form of equation (30) which on the line 6 =0 can be 
shown to be as follows: 


5 
~Mn(c? ie c’?) D) 1- a 
tou" “Te \afe\t aeee ae n 
(") (5) v\ =) (2) aa 
5° 
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Thus for the case of b=c = = 6 hence c’ =\/62+ c2- 6’2 oe 
4 


‘ b 
and at the point @ = 0, r= Me we have from equation (40) 
-Mn(c? — c’?) 2 
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1 0.22 
1-0.2 et 


Soa |? 
v( )? + 0.22 ( ra 5) ee 7 17 an) 
4 16 


=A Cee 


0.97017 - 0. : 
7) ae 98477] 


M > at 2) 
1.298 - MAGE Se 

= 
Hence equations (41) and (39) show that the central dipole, the 
central quadrupole, and the central octopole contribute 77.04%, 
12.41%, and 0.94% respectively of the total potential. Comparisons 

b 

of the undamaged double-layer of r= oe = with the total po- 
tential expressed in equation (36) show that the damage of the 
apex portion increases the importance of the central dipole and de- 
creases the importance of the central quadrupole and the central 
octopole. 

Similarly, any axi-symmetrically damaged double-layer can be 
analyzed by such superposition process. 

Disturbance potentials due to the presence of a spherical bound- 
ary. So far we have considered the potentials ¢, due to various 
current generators in a conducting medium of infinite extent. If a 

C) 
spherical boundary r=a>/ with the boundary condition “fo 


= =0 is imposed around the current generator, the disturbance 
, 

potential ¢, can be expressed in terms of the multipole strengths 
as follows (see Yeh, Martinek, and de Beaumont; 1958): 


$, =D, o + D, bf? + D, gf? + 
Cag Qe ee +8. or + 
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where 
2y 2 
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Thus for the eccentric radial dipole of strength pw, at (0, 0, 5) 
the disturbance potential is 


22(527 —3r?) 
2 
1 


fhe of BY LS 3/b\/7r\2 (44) 
“a ‘\7~]} |2 cos 6+ —|—]\—} (8 cos? @— 1) + 
r a 2\r/\a 
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is ah a 
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By using the interior sphere theorem of Ludford, Martinek, and 
Yeh (1955) the exact disturbance potential corresponding to ¢),, 


in equation (14) can be expressed as follows (see Yeh, Martinek, 
and Ludford, 1955): 


2 
a 
—a*n, (. _ -) 
a 


2 (45) 
ies \8 |/a2\8 cos 6- a 2 
r? a br a® a?\2)+ ” Or 
—2\—Jcos 6+ | — 
br 
At the point 0=0, r=a= at on the boundary, we have 
hi 
= Pal 1-(55) 1 
din = —et (1)* {| —— ] >——*G + Gl = 
fog r? (1) (0.4) ( a ) ei 0.4 
Si as 2S 
0.4 (46) 


mn 
+ 4.4444 2 =o 1260 Pou, ? 


and equations (44) and (46) show that the central dipole, central 
quadrupole, and central octopole contribute 45%, 27%, and 14.47% 
respectively of the total disturbance potential. At this same point, 
equations (13), (14), (44), and (46) together show that the central 
dipole, central quadrupole, and central octopole contribute 41.54%, 
27.69%, and 15.51% respectively of the total potential Pou, oa ips 
These check with the results of Frank (1953, Fig. 4), namely, that 
the total potential on the spherical boundary along the axis due to 
a radial eccentric dipole at b = 0.44 is 2.407 times the potential 
due to the central dipole. In addition, the computations here indi- 
cate that the potentials due to the central quadrupole and central 
octopole are 0.667 and 0.373 times respectively that due to the 
central dipole; and the combination of the higher-order central 
multipoles contributes a potential equal to 0.367 times that of the 


central dipole. 
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b 
At the point 6=0, r=a= aoe on the boundary, we have from 


equation (45) 


—~pE 1 A Hy 
aed (|| ee 47 
Gti peer? (dla aT) Z | : Pou, ? Ce 


and equations (44) and (47) show that the central dipole, central 
quadrupole, and central octopole contribute 71.11%, 21.33%, and 
5.69% respectively of the total disturbance potential. At this same 
point equations (13), (14), (44), and (47) together show that the 
central dipole, central quadrupole, and central octopole contribute 
68.57%, 22.86%, and 6.40% respectively of the total potential 
Pou, + Pip. And at the origin r = 0 the central dipole alone can 
account for the total disturbance potential due to the eccentric 
dipole in a spherical conductor. 

For the eccentric double-layer with an axi-symmetric circular 
rim of radius c in the plane z = 3, the disturbance potential is 


a 


Q2 3(327 — 1?) 
Zz Petes & 5 SNS Sees 
diy = 7Mc 73 + 2nMc* b rat 
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By the use of the interior sphere theorem the exact disturbance 


potential corresponding to ¢,y in equation (30) can be expressed 
as follows (see Yeh, Martinek, and Ludford, 1955): 


Pim = 2M — K(k) + 


my 


(b? + ¢2)7 E cos 6 — a2b/(b? +c?) 
a 
(49) 
A(e, F) — Oy, 


veaee ®,y is the value of ¢,y + ®,y at 7 = +0 and is a constant 


7,2 =17? + 2ra?(c sin 6-5 cos 0)/(b? + c*) + a4/(b? + ©), 


ECCENTRIC DIPOLE AND DOUBLE-LAYER 55 


FP = r2 —~2Qra?(c sin 6+ b cos 6)/(b? + e7\+ a*/(b2 4 07), 
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eee = 
og eS 
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e=F(y, k’) (51) 
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and 
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For the case of 6=c and at the point 6=0, r=a= ae on the 
boundary, we have 
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and equations (54) and (48) show that the central dipole, central 
quadrupole, and central octopole contribute 55.36%, 33.22%, and 
13.29% respectively of the total disturbance potential. At this 
same point equations (35), (29), (54), and (48) together show that 
the central dipole, central quadrupole and central octopole con- 
tribute 52.52%, 35.02%, and 14.71% respectively of the total po- 
tential doy +¢,y- This shows that the total potential on the 
spherical boundary along the axis of symmetry is 1.904 times the 
potential due to the central dipole. This checks with Frank’s 
computations (1953, Fig. 4) nicely, although he did not have a 
closed expression for the total potential. In addition, the compu- 
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tations here indicate that the potentials due to the central quadru- 
pole and the central octopole are 0.667 and 0.280 times respec- 
tively that due to the central dipole; and the total potential due to 
all higher-order central multipoles is equal to —0.043 times that of 
the central dipole. 


For the case of 6 =c and at the point 0=0, r=a= ‘o on the 


boundary we have 
‘A Mro2 2 ~ [(0.2)? + (0.2)2]7 
he p2 (0.2)? 1 is 
; 0.2 
~ (0.2)? + (0.2)? 


eed, ac?) ©) 


Mne2 2 ~~ [(0.2)2 + (0.2)2]7 — 0.2 
2 (0.2)? 1 e 


Mnc? 


2.72391 


anes 1.825 doy » 
and equations (55) and (48) show that the central dipole, central 
quadrupole, and central octopole contribute 73.43%, 22.03%, and 
4.41% respectively of the total disturbance potential. At this same 
point equations (36), (29), (55), and (48) together show that the 
central dipole, central quadrupole, and central octopole contribute 
71.15%, 23.72%, and 4.98% respectively of the total potential 
Pom + Pim: 


Cc 
For the case of 6 = 0 and at the point 0=0, r=a= ie on the 


boundary we have 


3 Mac? 20.2 1 
—— =a Serna Pe Ta va Pee he aa 
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and equations (56) and (48) show that the central dipole, central 
quadrupole and central octopole contribute 101.98%, 0%, and 


(56) 
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~2.04% respectively of the total disturbance potential. At this 
same point equations (37), (29), (56), and (48) together show that 
the central dipole, central quadrupole, and central octopole con- 
tribute 102.31%, 0%, and -2.39% respectively of the total potential 


dom + Pim: : 
"The disturbance potentials and total potentials on the spherical 


boundary due to a partially damaged double-layer can be evaluated 


in a similar manner. 


Conclusions. It has been shown in a previous paper (Yeh, 
Martinek, and de Beaumont, 1958) that any current generator in a 
volume conductor can be represented by a combination of central 
multipoles. whose strengths can be determined by a sufficient num- 
ber of surface potential measurements. In this paper the multipole 
representations of an eccentric dipole, an eccentric axi-symmetric 
double-layer, and a partially damaged axi-symmetric double-layer 
are obtained. “The results show that within the range of parameters 
corresponding to human subjects the quadrupole and octopole can 
contribute significantly to the surface potentials compared with 
the dipole. 

For both the eccentric dipole and the double-layer, the strength 
of the central dipole indicates the resultant strength of the genera- 
tor only, whereas the strengths of higher-order multipoles can re- 
veal other refined aspects of the current generator. For instance, 
in the case of an eccentric radial dipole equations (10) and (11) 
show that the ratio of the quadrupole strength Q_, =2 by, and the 
dipole strength D, = », is equal to twice the eccentricity } of the 
eccentric dipole. Hence we have 


(57) 


and the current generator as an eccuntric radial dipole is thus 
completely defined. The same idea also applies to non-radial ec- 
centric dipoles. A paper dealing with the eccentricity determina- 
tion of a general dipole in a volume conductor by means of the 
multipole strengths will be published shortly. 

In the case of an eccentric double-layer with a circular rim, the 
ratio of the quadrupole strength Q|. =27Mcb and the dipole 
strength D, = Mc? is again equal to twice the eccentricity 5 of 
the double-layer as indicated by equation (57). Furthermore, the 
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ratio of the octopole strength O,,, = nite? (A b? ~c”) and the 


dipole strength D, = 7Mc? is equal to 3b? - =o, which renders 


the radius of the circular rim c as follows: 


D s\o)- (58) 


With the radius of the rim thus determined the uniform dipole mo- 
ment density M can then be obtained as 


ae (59) 


and thus the current generator as an axi-symmetric double-layer is 


completely defined. 
During the process of the evaluations it has been revealed that 


in the multipole representation method used, there are degrees of 
freedom of the current sources which do not give rise to any po- 
tential outside the domain of singularities. Thus a double-layer on 
a part of a sphere and a circular double-layer disc admit different 
_ multipole representations giving exactly the same potential distri- 
bution outside the domain of singularities. A non-redundant though 
less graphic multipole representation is possible by employing the 
_ spherical harmonic expansion for potentials. Formulation of such 
a representation will be carried out next. Also under study is the 
formulation of an eccentric multipole representation which com- 
_ bined with the determination of the eccentricity of a general dipole 
_ mentioned above should give a much sharper picture of the current 
generator in a volume conductor than anything available at the 
present time. 


_ The authors wish to thank Dr. J. Watt, Dr. K. S. Cole, Dr. R. Pp 
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research, and Mrs. S. F. Kandel for editorial assistance in the 
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The present-day practices of electrocardiography and vectorcardiog- 
raphy are based upon the theory that the surface potential differences can 
be assumed to be due to a single dipole inside the body. It is shown in 
this paper that a dipole cannot account for all the surface potentials due 
to realistic current generators, and hence the determination of the current 
generator from surface potential measurements based upon such a theory 
will lead to inconsistent representations of the heart for one and the 
same subject. 

To demonstrate this point two eccentric dipoles of different strengths 
and locations representing two muscle fibers are taken to be the current 
generator in a homogeneous spherical conductor. The exact surface 
potentials are then expressed by means of the ‘‘interior sphere theorem’’ 
of the authors. With these expressions the magnitude, direction, and 
location of the resultant dipole are determined by the method of D. Gabor 
and C. V. Nelson (J. App. Physics, 25, 413-16, 1954). The surface 
potentials due to this resultant dipole are again exactly expressed by 
means of the ‘‘interior sphere theorem’’ and compared with those due to 
the eccentric dipoles assumed. It can be seen that the differences can 
be considerable. 

It is suggested that the multipole model of the authors (Bull. Math. 
Biophysics, 20, 203-16, 1958) be used as a more accurate and the only 
unique representation of the heart. 


Surface potential due to two eccentric radial dipoles in a homo- 
geneous spherical conductor. Consider a radial dipole p,= 
ip sina +kp cos a at the point a, =ia sin &+k@ COS a. The 


*This investigation was supported by the National Heart Institute 
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electrical potential due to this dipole in an infinite homogeneous 
medium of conductivity y is (see, for example, Yeh, Martinek, and 


Ludford; 1955) 


uw sin O(a —@ sin &) +p cos &(2-@ Cos O) 


V, (1, 8,9) = ne to 
4ryl(e -w sin a4)? + y7+(z2—@ cos a)*]? 


rn r(sin & cos Y sin 6+ cos & cos #)-—a@ 


» (1) 


a bed 


47y [724 2 —-2ra (sin 4 cos y sin 0 + cos 4 cos @)] 


where r, 0, are spherical polar coordinates which can be related 
to the rectangular coordinates 2, y, 2 as follows: 


a=rsindcos 9, y=rsin@sin 9g, 2=rcos 6. (2) 


If a spherical boundary r= @> @ is imposed around this dipole, 
the disturbance potential Vite 6, ~) can be obtained from the in- 
terior sphere theorem of Ludford, Martinek, and Yeh (1955) and has 
been derived in Yeh, Martinek, and Ludford (1955).: In this case 
it reads 
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. On the spherical surface the total potential distribution V, , (a, 6, 9) 
is then from equations (1) and (3) as follows: 


Ve (a, 6, ~) = vo (a, 0, ?) zr yy (a, 0, P) = 


n(a? — w ”) 


3 
4nryola?+o” —2a(sin 4 cos 9 sin + cos 4 cos 6)? (4) 


ph 
4nyaw- 


Consider now another radial dipole», =i sin (—-%) + kp cos (~ a) 
at the point w, =i sin (—&)+ko cos (—a). The total potential 
distribution V_,(a, 6, p) can be obtained from equation (4) with « re- 
placed by —a as 


V_o(4; 0, ~) = 
n(a? — w”) 
30 
47yala? + w* -2a@(- sin 4 cos P sin 9+ cos a cos 6)}? (5) 
se aes = 
4nyao 


When these two radial dipoles exist simultaneously in the spher- 
ical conductor the total surface gotential distribution will be 


V (a, 9, 2) = Vy gas 9; ~) + Via, 9, p). (6) 


Determination of the ‘‘equivalent” dipole. The present-day 
practices of electrocardiography and vectorcardiography imply that 
the surface potentials due to the heart as a current generator can 
be assumed to be due to a single dipole inside the body. With 
this in mind, we shall try here to determine the ‘‘equivalent” 
single dipole producing the surface potentials V(a,0,) due to 
the two eccentric radial dipoles considered above. One way for 
such a determination has been published by Gabor and Nelson 
(1954) based upon the concept of resultant dipole from surface 
measurements. 

Suppose the resultant dipole M =iM, +jM, +kM, is at the 
point P =iX+jY¥+kZ. For the case under consideration it is 
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obvious that we should have 


M,=M, =9, 
M, = 2p cos G, (7) 
X=Y=0. 


It now remains to determine the location Z along the z-axis. The 
method of Gabor and Nelson (loc. cit.) involves the evaluation of 
two surface integrals in the expression 


ya 277 7 ; 
z-— 1s | ag { V(a, 6, p) cos? 6 sin 0d6 - 
2p Cos a Jy 0 
ya® 277 7 
| cos? pd { V (a, 6, p) sin® 6d8. (8) 
2p COS O J 0 


Such evaluation becomes prohibitive, if not impossible, even for 
the simple case considered. However, if the dipole is indeed 
‘fequivalent,’’ it should produce the same potential distribution as 
V (a, 6, ~) in equation (6) on the surface. Hence we need only to 
match the potentials on one of the surface points in order to deter- 
mine Z. 

The surface potential distribution due to this ‘‘equivalent’’ 
dipole can be obtained by putting a=0, w=Z, and p=M,= 
2cos & in equation (4) as 


Wea 0: 9) oe eats a ie 
4nyZ(a2+Z?~—2aZ cos 6)* 
2p cos & 
4nyaZ : 
Let us now match the potentials at the point r= a, 6=0, that is 


(ie (a, 0, p) = V(a, 0, ~) = Vio (as 0, p) + V_i(a, 0, 9). (10) 


(9) 


Substituting equations (4), (5), and (9) into equation (10), we have 
cos a(a?-Z*) cos a 
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Dertieetcue =, 
( ~ =] cos 1 + i -~ 2 (=) cos ae 
a a a 


geet 4 
which is a quadratic equation in — and can be solved to express 
a 


Z in terms of a, w, and a. 
Consider, for instance, the case of w=0.4a@ and a = 30°. 
Equation (12) becomes 


: Z 
ag 2.5 — 0.4 2.5 
= - — = 4.70709, (13) 
At eee ‘3 tz ay 
( -=| B fi +0.18-08 aa aks 
a 2 ) 9 
or 
Z\2 Z 
est Sios 4 —} 4086266 =O, (14) 
a a 
Hence 
F 
—= 0.89378 +,/(0.89378)? — 0.36266 = 0.89378 + 
a (15) 
1.55422, 
0.66044 = ea 


from which we take Z = 0.233344 as the solution to our problem. 

Similarly, if we match the potentials at the point r= a, 0 = 180°, 
we obtain Z = 0.31149a@ instead. This already shows that an 
‘“equivalent’? dipole giving the same potential at one point 6 =0 
cannot account for the same potential at the other point 0 = 180°. 


Surface potentials due to the determined ‘‘equivalent” dipole. 
The surface potential distribution due to the determined **equiva- 
lent’? dipoleM =k 2p cos 30° =kV/3patthe pointP =k- (0.23334 a) 
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can then be obtained from equation (9) 


v3 p 
Vega 6, ~) gi Geer a eo 
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1 — (0.23334)? 
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{1 +. (0.23334)? ~ 2 (0.23334) cos 6}? 
Mm 0.55853 
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From equations (4), (5), and (6) the surface potential distribution 
due to the two eccentric dipoles p, =ip sin 30° + kp cos 30° = 
i 0.5» +k 0.86603 p at w, =i 0.4 asin 30° +k 0.4 acos 30° = 
i0.2¢+k0.34641a and p,=ip sin(-30°)+kp cos (-30°) = 
-i 0.54 +k 0.86603 at w,=i(0.4a)sin(-30°)+k (0.42) cos (-30°) = 
~i10.2a@+k0.34641 a is 


Sod eee 
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{1 + (0.4)? — 2(0.4)(-0.5 cos 9 sin 6 + 0.86603 cos ay} | (17) 
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~ 0.089 . 


(1.16 + 0.4 cos p sin @ — 0.69282 cos at 
_ It is obvious that equations (16) and (17) do not render identical 
potentials everywhere on the spherical surface. Equation (16) is 
a function of @ only and is constant for all y’s at each 6, whereas 
equation (17) is a function of both 6 and —, being symmetric with 
respect to p = 0 and ¢ = 180° as well as with respect to p = 90° 
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and ~ = 270°. Therefore it is equally obvious that no single di- 
pole can accurately reproduce the surface potentials due to two 
dipoles. Thus we arrive at the same conclusion as Okada (1958) 
who also presented graphically in his Fig. 2 the discrepancies 
between the surface potentials indicated in equations (16) and 
(17). Incidentally, his Fig. 2 seems to show that Vz g (a 6, p) and 
V (a, 6,9) are identical both at 6=0 and 6= 180°. But our re- 
sults show that if they are identical at 6 = 0 there is a difference 
of approximately 10 per cent at 6 = 180°. 


Determination of the central dipole from surface potential dif- 
ferences. We can demonstrate the inconsistency of the single 
dipole concept in another way. In the theory of the conventional 
lead systems in electrocardiography, any three sets of surface 
potential differences can determine a central dipole representing 
the current generator. Take, for example, Wilson’s tetrahedron 
system. With the consideration of the ‘‘central terminal’’ we may 
determine the central dipole strength by potential values at three 

- -a -—2,2 
points, say P,, eee = = ee 0, = ay , and 
P,, (0,4, 0).: For the numerical example considered, we have from 


equation (17) 
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Then the components of the central dipole strength are (see, for 
example, Yeh, Martinek, and de Beaumont, 1958): 
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Thus we obtain a central dipole representation which is not sym- 
metric, either with respect to the xz plane or with respect to the 
y 2 plane. 

Suppose now we realize the fact that the current generator is 
symmetric with respect to the zz plane. Then we must have 
D =0 for the central dipole representation. We next proceed to 
determine D, and D_ from potential differences among three points 
in the x2 plane corresponding to the apices of an Einthoven tri- 
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surface potentials at these points are, from equation (17), 
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Then the D, and D, components of the central dipole strength are 
(see, for example, Yeh, Martinek, and de Beaumont, 1958) 
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Thus we obtain a central dipole representation which is different 
from that indicated by equation (19) and which is symmetric with 
respect to the zz plane as well as the yz plane, as it should be. 
But this is by no means the unique representation. 

In general, unless the surface potentials are actually due to a 
central dipole, the measurements on different sets of points on 
one and the same subject will lead to different central dipole 
representations of the current generator, since no single dipole 
can account for all the surface potentials which, for well-known 
reasons, are simply not due to a single dipole. 


Conclusions. The heart as a current generator in a homogene- 
ous conductor can be considered approximately as a combination 
of a number of discrete and distributed dipole elements in a po- 
tential field governed by the Laplace equation. In this paper it is 
shown that even for the simplest type of dipole combination, i.e., 
that of two eccentric radial dipoles, a single dipole cannot ac- 
count for all the potentials on the surface of the conductor. Hence 
the present-day mathematical-physical methods in electrocardiog- 
raphy and vectorcardiography which are based upon the single 
dipole concept will lead to inconsistent representations of the 
heart for one and the same subject. 

Even if one would insist on employing the single dipole con- 
cept, there are satisfactory means to determine only the strength 
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of the ‘‘equivalent’’ dipole, but not the location. The three- 
dimensional expressions of Gabor and Nelson (loc. cit.) for re- 
sultant dipole strength and location are formally simple but can 
rarely be carried out in practice. 

It is suggested that the multipole model of Yeh, Martinek, and 
de Beaumont (1958) be further explored, as it is a more accurate 
and the only unique representation of the heart if merely measure- 
ments on the body surface are given. In that model, a sufficient 
number of central multipoles can account for all the measurable 
surface potentials. The strengths of these central multipoles can 
then be uniquely determined from a sufficient number of surface 
measurements, regardless of where these measurements are taken. 
Hence a consistent representation of the heart on one and the 
same subject is assured. 
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The DNA-protein coding problem is given a general algebraic formu- 
lation, the Gonsequences of which are then explored by standard mathe- 
matical methods. To keep the treatment self-contained, the mathematical 
techniques to be used are explained in detail. It is demonstrated that 
there exist @ priori a countably infinite number of different abstract 
DNA-protein codes, thereby showing that inductive attempts to construct 
such a code will most likely be fruitless. A notion of ergodicity is then 
introduced, which imposes a number of restrictions on the admissible 
codes, and, in fact, these considerations enable us to derive a small 
portion of a code which, if our hypothesis of ergodicity is correct, must 
occur in nature. Finally, we discuss briefly the problem as to whether 
there can exist more than one DNA-protein code in nature. 


The purpose of this paper is to discuss some of the quantitative 
aspects of the mechanism of the coding of protein molecules by 
the DNA of a living cell. Previous studies (Gamow et HL LoD), 
1956) have mainly emphasized the statistical aspects of the cod- 
ing problem; our approach, on the other hand, will be of a purely 
algebraic nature. In particular, we shall be concerned primarily 
with the following three problems: (1) How many different codes 
between protein and DNA are a prior possible? (2) Is the statis- 
tical character of the code of the particularly simple type which 
is usually defined as ergodic? (3) Can there exist more than one 
DNA-protein code in nature? Since the algebraic techniques that 
we shall employ are somewhat unfamiliar to non-mathematicians, 
we shall provide a complete discussion of them as they are 
required. 

We shall, as usual, consider that an arbitrary DNA molecule is 
a polymer of four different elementary substances (nucleotides), 
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while a polypeptide chain is regarded as a polymer of 20 different 
elementary substances (amino acids). According to the well- 
known Watson-Crick hypothesis (1953), the DNA molecule is 
thought to consist of two helically intertwined complementary 
strands, on each of which the constituent nucleotides are arranged 
in a sequence, or linear array. We may restrict our attention to 
these individual strands (since a single such strand determines 
the entire molecule), or else we may order the nucleotides of a 
given DNA molecule into a sequence by listing them serially as 
we actually traverse the molecule, as in the scheme of Gamow 
(1954); in any event, the polynucleotides will occur in a definite 
sequence along the length of the DNA molecule. It is generally 
supposed that the sequence of nucleotides in a molecule of DNA 
will uniquely determine the sequence of amino acids in a polypep- 
tide chain; the DNA-protein coding problem is, then, concerned 
with determining the precise manner in which sequences of nucleo- 
tides can determine sequences of amino acids. While this may be 
thought to be largely a biochemical question, it will be seen that 
the purely mathematical discussion provides, at least, some in- 
teresting insights. 

We shall now introduce the basic algebraic tools which will be 
required in the sequel. Let us suppose that S ={0,,%5,+++,0,} 
is a finite set, and let us denote by M the totality of all finite 
sequences of elements of S; we remark explicitly that we allow 
repetition of any element of S. Each such sequence will be de- 
noted by juxtaposing its terms; thus, an arbitrary element of M 
will be of the form 


U= eat tes ee Oi 
where a, €S for 1S’ <n An element of M will be called a word; 
the number of elements of S (distinct or not) in the word will be 
called the length of the word. Thus, the length of the word dis- 
played above is r, even if all the a; , are the same. 
We may define a law of etmapdaltibe in M as follows: let w, = 

Fey Oi 030, wy =e 2 -».Q. be arbitrary words of M. We 
define the phdiee of w, and We as the word 


Rg Uh = Q. Oe atte ees Qa. » @ ee eC} 
tiers ty ty tr te nF 


ye 


DNA-PROTEIN CODING (6) 


It is readily verified that this composition operation satisfies the 
associative law; i.e., given any three words w,, w,; W, of M, we 
always have (w, w,)w, = w,(w, w,)-' It is not commutative; i.e., 
W, Wy # Wy Uy) in general. If we introduce the notion of the word 
of length zero (called the empty word) which we shall denoto by 
e, then we see that for an arbitrary word we M, we have 
ew = we =U, 

so that e is a neutral element for our law of composition. It is 
clearly the only neutral element, for if e’ were another, we should 
have e’= e’e=e. 

Thus, we see that our law of composition has converted the set 
M into an algebraic system, satisfying the following conditions: 

1. M is closed under the composition; i.e., wy, W,¢M implies 
that w, w eM. 

2. The composition is associative. 

8. There exists a neutral element, e. 

An algebraic system which satisfies these three properties is 
called a monoid (see e.g., Bourbaki, 1951, ‘Chevalley, 1956). The 
monoid M which we have constructed above is termed the free 
monoid generated by the elements of the set S; these elements are 
called the free generators of M. 

We shall show that these free monoids are the proper algebraic 
vehicle for the study of those aspects of the DNA-protein coding 
problem in which we are interested. The motivation and biological 
significance of the word ‘free’? in the definition of M, will be 
explained below, following the development of those aspects of 
the theory of monoids which will be necessary for the relevant 
arguments. 

The first important notion is that of a submonoid of a monoid. 
If M is a monoid, then a subset M’CM is called a submonoid of M 
if it satisfies the following two conditions: . 

Si: If w, w’eM4, then wu’eM’ (i.e., M’ is closed under the com- 
position in ™).' 

S9: If e is the neutral element of M, then ee M’. 

Explicitly, (S2) requires that a submonoid of M have the same 
neutral element as M itself. That this is, in general, a real re- 
striction, can be seen from the following example: take M to be 
the set of integers reduced modulo 12, with the composition opera- 
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tion taken to be multiplication mod 12; let M’ be the set {0, 3, 6, 
9}. It may be verified that M’ is closed under multiplication mod 
12, and that 9 acts as a neutral element for M’. Hence, M’ is a 
monoid; it is not, however, a submonoid of M, since the neutral 
element for M is 1. 

The next concept which must be introduced is that of homo- 
morphisms of monoids. If A and B are monoids, then a mapping f 
from A to B (which we shall henceforth symbolize by f: A > B) is 
called a homomorphism if it satisfies the properties: 

H1: f(a, @,) = f(a,) f(a.) for every a,, a, A. 

H2: If e , denotes the neutral element for A, and e, denotes the 
neutral element for B, then f(e,) = ep: 

Let us illustrate this concept by constructing several concrete 
-examples of homomorphisms, which will be of considerable im- 
portance in the sequel. Let us denote by Z* the set of all non- 
negative integers; i.e., Zs (0; 8 Se eis readily verified 
that the usual operation of addition of integers converts Z* into a 
commutative monoid, the neutral element of which is 0.' Let M be 
the free monoid generated by the set S={d,,++-, of. ' We have 
already associated with each word weM a non-negative integer, 
its length, which we shall denote by A(w). It is readily verified 
that this mapping \:M — Z* satisfies the properties (H1), (H2), 
and is hence a homomorphism between M and Z*.. We shall refer 
to this homomorphism as the length mapping of M. In a similar 
fashion we can define further mappings \,:¢@ — Z*, where A; (wv) = 
the number of times the generator 4, appears in the word w. Once 
again it is easy to verify that for each 7, 1<i<n, the mapping 
A, is a monoid homomorphism between M and Z*; the mapping d; 
will be referred to as the ith partial length mapping of M. The 
partial length mappings are seen to satisfy the property 


> Ag (to) = A(w) 
iml 
for each word weM. . 
As a final property of homomorphisms, we may remark that if 
f: A — B is a homomorphism of monoids, then the image of A 
under f (i.e., the set of all elements b¢B such that $= f(a) for 


some ae A) is a submonoid of B; this follows readily from the 
definitions. 
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Like all other mappings, a homomorphism f: A — 5B may have 
the property that it is onto, in which case it is called an epi- 
morphism; it may be 1-1, in which case it is called a monomor- 
phism; or it may be both 1-1 and onto, and is then called an iso- 
morphism. Two isomorphic monoids are abstractly identical. It is 
easy to check that the length mappings we have defined above are 
epimorphisms but not monomorphisms. Much of the discussion 
below will involve the construction of monomorphisms between 
free monoids; it will in fact be seen below that the concept of 
monomorphism is precisely the one required for our formulation of 
the coding problem. 

' The statement that the mapping f: A — B is a monomorphism is 
equivalent to the statement that A is isomorphic to a submonoid of 
B. Under these circumstances, we shall sometimes say that A can 
be embedded (or isomorphically embedded) as a submonoid of B. 

The next notion which will be introduced is that of a quotient 
monoid; we now enter into a brief discussion of this important 
concept. A more detailed discussion of quotient structures may 
be found, e.g., in Chevalley (1956). 

It is well known that an equivalence relation defined on a set 
T will partition T into a family of disjoint, non-empty subsets, 
which are called equivalence classes. If we denote the given 
equivalence relation by F, then the set of equivalence classes of 
T under R is called the quotient set of T under R, and is denoted 
by T/R. The mapping 7:T — T/R, which associates to each 
element te 7 its equivalence class, is called the natural mapping 
of T onto T/R. 

Now let 7 be a monoid M. Given an equivalence relation # on 
M, we can always use the above construction to obtain a quotient 
set M/R. However, we may ask what conditions must be imposed 
on R in order that M/F will also be converted into a monoid in a 
natural fashion. We denote the statement ‘7 is equivalent to 6 
under R” by akb, and shall prove that the required condition is 
the following: 


aRa’ and 6Rb’ implies (ad) R (a’d’) (1) 
for all elements a, a’, 5b, 6’ in M; an equivalence relation on M 


which satisfies this condition is called compatible with the multi- 
plication in M.. We assert: if R is any equivalence relation on a 
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monoid M compatible with the multiplication on M, then M/F is a 
monoid, and the natural mapping 7:M — M/F is an epimorphism. 

To show this, we must exhibit the law of composition on M/R 
which converts it into a monoid; that is, given two equivalence 
classes A, BeM/R, we must define a product class AB. To do 
this, let us select ae A, b¢B as arbitrary representatives of their 
equivalence classes. Now a and 6 are elements of M, which is 
already a monoid; hence the product ad is defined. We define the 
product AB of the equivalence classes of a and 6 as the equiva- 
lence class of the product a5. Conditions (1) ensure precisely 
that this multiplication in M/R is well defined (i.e., independent 
of the choice of representative elements from A and B). We leave 
it to the reader to verify that this composition is associative, that 
there is a neutral element E£, which is precisely the equivalence 
class of the neutral element e of M, and that the natural mapping 
7 is indeed an epimorphism. The monoid M/F so defined is called 
the quotient monoid of M modulo F. 

We may now explain the biological significance of the struc- 
tures we have introduced. Our point of departure is the observa- 
tion that we may consider the totality of possible nucleotide 
chains (i.e., the totality of possible DNA molecules) as possess- 
ing the algebraic structure of the free monoid on four free genera- 
tors (which correspond to the individual nucleotides), and the 
totality of polypeptide chains as possessing the structure of the 
free monoid on 20 generators (corresponding to the fundamental 
amino acids). The biological problem of associating to each 
polypeptide a DNA molecule which ‘‘codes’’ it then reduces to 
the problem of constructing monomorphisms of the free monoid on 
20 generators into the free monoid on four generators. (It should 
be remarked that the numbers four and 20 are in no way essential 
to our arguments; they are merely numbers chosen, in conformity 
with current literature, for the construction of suggestive examples 
illustrating our techniques). 

Let us pause for a moment to examine the mathematical and 
biological significance of the word ‘‘free,’’ as applied to the 
monoid we constructed on p. 73. A free monoid generated by a 
set S is, in a sense, the largest possible monoid which S can 
generate. The free monoid is, roughly speaking, characterized by 
‘the fact that there exist no relations in it of the form W, = We; 
unless these words are identical. In a monoid which is not free, 
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we may still be able to find a finite set S such that every element 
of the monoid is a finite sequence of elements of S, but there will 
exist relations of the above form between some of these sequences. 
In general, if we take a free monoid M and impose upon it a finite 
number of relations between the words of M, then it is not hard to 
show that, if we associate to each word weM the totality of all 
words «* such that w = w’, we obtain thereby a decomposition of M 
into equivalence classes, corresponding to a compatible equiva- 
lence relation R on M. The quotient monoid M/F so obtained cor- 
responds to a monoid generated by the same set S that generates 
M, but which is not free. In fact, free monoids can be charac- 
terized by the following universal property: every monoid is a 
quotient monoid of a suitable free monoid. Since we shall not need 
this result in the sequel, we shall not enter into its proof here. 

In biological terms, the significance of our algebraic structures 
is that they represent the currently accepted biochemical hy- 
potheses concerning the structure of nucleic acids and polypep- 
tides. Thus, for example, our representation of the totality of 
nucleic acids, and the totality of polypeptides, by free monoids 
embodies the biochemical assumption that all conceivable poly- 
mers of nucleotides and amino acids are capable of existing 
physically. Likewise, our assumption that the DNA-coding process 
is equivalent to the choice of monomorphisms between appropriate 
free monoids embodies a generalized one gene-one enzyme hy- 
pothesis, since once a definite code (i.e., monomorphism) is 
chosen, then, by definition, distinct polypeptides must correspond 
to distinct DNA molecules. The question, however, as to whether 
more than one such monomorphism can exist in nature will be 
taken up later. 

Further, our formulation of the coding problem in terms of 
monomorphisms between free monoids already restricts us to the 
consideration of the comma-less codes introduced by Crick, 
Griffith and Orgel (1957). To see this, let us introduce the fol- 
lowing definition: a word w of a free monoid M is said to be 
embedded in a word w of M if there exist words v,, v,eM such 
that W=v,wv,. If now M,, M, are the free monoids generated by 
sets of m and n elements respectively, and /:M, — M, is a 
monomorphism, then, by definition, Mis isomorphic to a sub- 
monoid of M,. In fact, if {b,, Borers 6} is the set of generators 
of M_, then M is isomorphic to the submonoid of M, generated 
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by the words f(5,),+++,/f(5,,).' Now suppose that w,, wv, are 
words of M, such that f(w,) is embedded in f(w,).' That is, 
there exist words v,, ¥,¢f(M,,) (not merely in M_) such that 
fGe5) = vf) ey But since v,, v,€f(M)s there exist, by 
definition, words u,, uv ,¢”@, such that v, = f(u,)> %_ = f(ug)- 
Therefore, we have f(w,) = f(u,)f(w,) f(u,) = fu, wv, u,)’ This 
argument shows that if w,, w, are words of M_ such that f(w,) is 
embedded in f(w,), then w, is already embedded in w,- But this 
is precisely the definition of comma-less code given by Crick 
et al (loc. cit.), adapted to our terminology. 

Therefore, in the light of what has been said above, we may 
observe that our formulation of the DNA-protein coding problem is 
not an abstraction from the generally accepted biochemical situa- 
tion, but a precise translation of this situation into algebraic 
terminology. The advantage of this formulation is that the power- 
ful methods of modern algebra may now be directly applied to the 
problem. As a result, all of the consequences of the mathematical 
theory are directly applicable to the physical situation, and.any 
discrepancy between the mathematical development and actual 
biochemical observation must (unless the mathematical arguments 
contain an error) derive from incorrect hypotheses on the bio- 
chemical level. Further, the purely algebraic approach must neces- 
sarily yield the same results under these circumstances as a more 
empirical approach derived from stereochemical, analytical and 
X-ray diffraction data. 

We are now ready to undertake the argument itself. We shall 
proceed in full generality, but we shall attempt to motivate the 
somewhat involved constructions by means of relevant biological 
examples. Our first result (Theorem 1 below) asserts that there 
exist a countably infinite number of a priori possible DNA-protein 
codes; i.e., there exist infinitely many different monomorphisms 
which embed the free monoid on 20 generators into the free monoid 
on four generators. To do this, we will require a number of pre- 
liminary constructions and lemmas. The plan of the proof is as 
follows: we shall introduce, analogous to the notion of free monoid, 
the notion of a free group (a group, algebraically, is a special 
case of a monoid, satisfying, in addition to the monoid properties, 
the additional condition that each element shall possess a multi- 
plicative inverse); these groups are already known to possess the 
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properties which we would like to demonstrate in an arbitrary free 
monoid. Our first lemma will show that a free group can be con- 
sidered as a quotient monoid of a suitable free monoid. In Lemma 
2, we write down an important property of free groups (namely, 
that every sub-group of a free group is itself a free group), which 
we will not be able to show for general monoids until much later 
(Theorem 3). In Lemma 3, we show that, given an arbitrary sub- 
group of a free group, there is a canonical manner of passing to a 
submonoid of the free monoid of which the free group is a quotient. 
Finally, Lemma 4 asserts that the free group on two generators 
already contains the free group on n generators (for every integer 
n 2 2) monomorphically embedded in it. Theorem 1 then combines 
these lemmas to give us our result. 

We now proceed to the details. Let S be a finite set, and let 
us denote by S’ a set which can be put into 1-1 correspondence 
with S Let j: S$ —, S’ denote a definite 1-1 correspondence be- 
tween S and S% Let ™M be the free monoid generated by the set 
S US. In order to obtain a group from M, we must have some 
means of constructing inverses; the plan is, roughly speaking, to 
consider j(%) to be « *. With this in mind, we observe that we 
can extend j to M by writing j(a8) =/7(8)j(a) for all a, BeS; 
this equation makes sense if we consider all the terms as ele- 
ments of M. This definition will lead to the usual property of 
inverses; namely, (48) 1=B 'a7'. 

Let us now introduce an equivalence relation into M as follows: 
we shall call two words of M equivalent if one of them can be 
obtained from the other by the removal of any number of factors of 
the form aj(a) or j(a%) a. Remembering that j(H) is to correspond 
to «71, then the equivalence of two words of M means that words 
which differ by the occurrence of the neutral element any number 
of times are not to be regarded as essentially distinct. It is 
easily verified that this equivalence relation, which we shall 
denote by R js is, in fact, compatible with the monoid composition 
in M.: Hence, we may form the quotient monoid M/R.. We now 
prove 

Lemma 1. M/Rf. is a group. 

Proof. It is only necessary to show that each element of M/R, 
possesses an inverse. The multiplication in M/R, is obtained, as 
explained above, by choosing arbitrary representative elements in 
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M from the equivalence classes of M/E, multiplying them in M, 
and taking the equivalence class of this product. Since this 
multiplication is independent of the choice of representative, we 
may choose this representative in a canonical way; we agree to 
choose the element of smallest length in each equivalence class. 
By what we have said above, this corresponds to choosing words 
which do not contain redundant combinations %j(a)(=aa~'). 
Then if A is a class in M/R,, and ae A is the element of smallest 
length in A, we take the class of j(a) and denote it by A=. CGE 
it then follows immediately from the definitions that AA~' = A~'A = 
E. Thus the lemma is proved. 

The group M/F, which we shall designate by G, is called the 
free group generated by the set S. It might be thought that this 
group is not uniquely defined, since it apparently depends upon 
the choice of a definite 1-1 correspondence j: S$ — S. However, 
it can be shown that choosing a different mapping j*: S ~ S’ will 
result in a group which is isomorphic to the one we have con- 
structed using j; the proof of this assertion is left to the reader. 
The construction we have employed above to obtain the group G 
appears as an exercise (p. 27) in Chevalley (loc. cit.). 

Lemma 2. Every subgroup of a free group is itself a free group, 
generated by some subset of the group. 

Proof. This result is known as the Nielsen-Schreier Theorem, 
and is a standard result of group theory. For a simple proof, see, 
e.g., Kurosh (1956), Vol. II, p. 33. 

We now prove the result which connects the theory of free 
groups with the theory of free monoids which has been developed 
above. 

Lemma 3. Let G be a free group, and let G be expressed in the 
form M/F, where M is a free monoid and R. is a suitable equiva- 
lence relation on M. If H is a subgroup of G, then there exists a 
free submonoid N of M such that H = N/R; 

Proof. From Lemma 2, we know that H is itself a free group, 
generated by a set 7, where 7 is a set of elements of G. That is, 
AeT implies that A is an equivalence class of elements of M. 
For each AeT7, let us choose the element ae A of smallest length. 
Let us consider the set consisting of all of these elements ay 
together with all elements of the form j(a) for each such a; this 
forms a subset of M. Since equivalence classes are disjoint, 
these elements must all be different; further, there is no element 
a in this set such that a = j(a), for this would imply that A? = E 


a a eh a OD 
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in G, whereas it is well known that no element of a free group, 
save the neutral element, can have finite order (cf. Kurosh, loc. 
cit.).' Hence, we may form the free monoid N on the set of these 
elements of M. But this is a submonoid of M, and it is further 
evident from the definitions of the concepts involved that H = N/K; 
as asserted. This proves the lemma. 

Our next result embodies the property of free groups which will 
be of interest to us. 

Lemma 4. The free group on a set consisting of two elements 
contains as a subgroup the free group generated by a countably 
infinite set of generators. A fortiori, the free group on two ele- 
ments contains as a subgroup the free group on a set of n elements 
for every positive integer n. 

Proof. This is also a standard result of group theory. We assert 
here that if G, is the free group on two generators, then the 
commutator subgroup (i.e., the subgroup generated by all elements 
of G, of the form zyx 'y', 2, yeG,) is free on a countably in- 
finite set of generators. Hence, for each integer n, we may choose 
any n generators of the commutator subgroup to exhibit a free 
group G_ on n generators embedded as a subgroup of G,.' For full 
details of the proof that the commutator subgroup actually pos- 
sesses this property, we refer the reader to Kurosh (loc. cit), 
especially Vol. Il., p. 36. 

Our first main result now follows without difficulty from the 
lemmas which we have proved: 

Theorem 1. The free monoid on 20 generators can be isomorphi- 
cally embedded as a submonoid of the free monoid on four genera- 
tors in a countable infinity of different ways. 

Proof. Let us agree to denote the free monoid on n generators as 
M,» We observe that, by virtue of Lemma 1, the free group G, can 
be regarded as a quotient monoid of the free monoid M,. Lemma 
4 tells us that the free group G, contains the free group G, as a 
subgroup, for each integer n- Using the construction of Lemma 3, 
we can find for each such free subgroup a free submonoid of M,. 
Hence, it follows immediately that we can, in particular, find 
countably many different ways of embedding M,, as a submonoid 
of M,. 

Theorem 1 thus establishes the fact that there are infinitely 
many DNA-protein codings which are possible 4 prior. Stated in 
full generality, Theorem 1 asserts that the free monoid M, con- 
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tains as a submonoid the free monoid M,, where n = 1, 2,..-, Np. 
We can use this result to obtain a still sharper statement: 

Theorem 2: Any free monoid M, (n= 1, 2,..., %)) can be em- 
bedded as a submonoid of the free monoid M,, 

Proof. By virtue of Theorem 1, it is clearly sufficient to show 
that M, can be embedded in M,.' The situation is now sufficiently 
simple for us to exhibit a specific embedding. To do this, sup- 
pose that M, is generated by the symbols a,, @,. ' Consider the 
words 


= 


We shall prove that the submonoid of M, which is generated by 
bis obse spel, is isomorphic to the free monoid M ye Let us con- 
sider the free monoid M, generated by the symbols «,, Xo Xs, 
a,, and let us define a mapping f:M, — M, by writing f(a,) = 5,, 
t= 1, 2, 3, 4. By the rule (H1) for homomorphisms of monoids, it 
is seen that it is enough to specify the behavior of the generators 
under the mapping f in order to determine the homomorphism. To 
show that f is a monomorphism, we must prove that if w,, w, are 
words of M, such that f(w,) = f(w,), then w, = w,. If 

Cos nor ple al oni 


= » ORS eee . 
W a ig O 5° 


2 s 


then, since fis a homomorphism, we have 
f(u ) = f(a; f(a; )- ae f(a; ) 
= b:. bg eee bi, 
and likewise, 
f (wa) = b; bi ee bys 
Hence, the equality of f(w,) and f(w,) means that 


b. bred, = 8; b. os Oe ° 


oy Jy Ja Js 


However, a glance at the definition of the b's, together with the 
fact that we are dealing with free monoids, shows that the only 
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way this situation can occur is if r=s and 2, =), for k= 1, 
2,+++,% . This proves the theorem. 

We may remark that Theorem 2 is an exact statement of the 
well-known fact of general coding theory, that any code whatsoever 
may be reduced to a binary code (i.e., a code utilizing just two 
symbols). This result is not of direct biological significance in 
itself; however, it will be needed to prove our next result, which, 
shortly, will be seen to provide us with interesting biological 
consequences. 

Theorem 3: Every submonoid of a free monoid (on a finite or 
enumerable set of generators) is itself free on a finite or enumera- 
ble set of generators. 

Proof. By Theorem 2, it suffices to show that every submonoid 
of M, is free on a finite or enumerable set of generators. 

Let Z* denote the set of non-negative integers (which is a 
monoid under the usual addition operation) and consider the map- 
ping A:M, > Z*, which associates to each element of M, its 
length. If N is a submonoid of M,, then the set A(N) is a set of 
non-negative integers. Since only the empty word of M, is mapped 
onto 0 by A, it follows that the set A(N) has a smallest non-zero 
element n,. The set [A~'(n,) ON] must then be a subset of the 
set of generators of N, since each element of this set is inde- 
composable (in N) into a product of non-empty words of smaller 
length. In fact, any element xeN such that we) < 2n, must 
also be a generator of N, for a similar reason. Let us designate 
by S, the set 


fa:aeN, r71(z) < Qn}. 


If S, generates N, the proof is finished, since S, must, in fact, be 
a Baits set (there can only be a finite number of words of a fixed 
length in M,).' Let N, denote the submonoid of N generated by 
83 clearly, N, is a on the elements of S,.' We shall suppose 
that N, #N. ee us write N’=N — Ng 3 No i is not a monoid in 
poser: “put the mapping d is still defined pr it. Letn, = min A(N’). 
It follows from our construction that n, 2 2n), and kee that 
every zeN’ such that A(z) <n, +n, is a generator of N. For if 
such an element x were not a generator, then 2 would be ex- 
pressible in the form z= 2, 2,, where a, a,eN. There are now 
_three possibilities: 
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1. @,) #, €S,.' But then aeN, , contradicting the hypothesis that 
1 


weN’=N-N, 

Qe @, wy eN. But obviously A(z) = A(@,) + A(z,).' But A(z, \ 
A(x ene and hence, A(z) 2 2n,.' However, 2n, >, +, con- 
tradicting the hypothesis that d(2) pod een 

8. '2,€S,, #,¢N% In this case, we pave i (2) >n, + 2n, (by the 
additivity of A), thereby again contradicting the er on the 
length of a. 

Hence, let 8, ={a:aeN%4 A(z) <n, +n,}' Let T, =8,U S,, and 
let N,. be the submonoid of N Ped oentea by 77. Tf N, =N, then 
we are done, since 7, is at most a finite Sst and bfodity N; 15 
free. If 7, does ot: generate N, we consider the set N’’= “NV - 
Nr, , and let n, denote min A(N**). : Here again, we pick out those 
words of N” onic by virtue of their length and the minimal prop- 
erty of the n,’s, 7 = 0, 1, 2, must be generators of N. We call this 


set S, and pereidee the set 7, = 7, US,. If the submonoid NV 
AEA N, we are done; sib onaine: ia can form the set Nea 2 
N —N*”.and proceed as before. It is easily seen that this process 
can continue for at most an enumerable number of steps, at each 
of which we determine a finite number of generators of N. . Hence, 
the resultant set T = 3 will be at most enumerable. Further, 
it is evident from our construction that every generator of N must 
fall into some set S,, and therefore, we have exactly N, =N. 
Hence, our theorem is proved. 

It may be remarked that Theorem 3 is the monoid analog of the 
Nielsen-Schreier Theorem for groups, which was quoted earlier as 
Lemma 2 above. It will be applied in our discussion of ergodicity 
in relation to the DNA-protein coding problem. Accordingly, we 
shall now undertake a brief description of general ergodic codes, 
as they appear in the light of the algebraic framework we have 
constructed above. 

In ordinary coding theory, to say that a code source has been 
specified is to say that two items have been given: (i) the set of 
symbols (alphabet) of which the messages emitted by the code 
consist, and (ii) the statistical structure of the source. This 
second requirement is fulfilled by giving some rule which governs 
the statistical behavior of the ensemble of messages emitted by 
the source. By far the simplest and most important such rule is 


WE Pa te % 


my Ey Ye 
At 


a ara A it ah ak 5 
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the following: every sequence of the source which is of suf- 
ficiently great length has the same proportion of letters of the 
alphabet used by the source. This property is called ergodicity 
(Khinchin, 1957). 

We shall now attempt to formulate this concept in algebraic 
terms. Let 1O,, Agyeeey a} be the alphabet of an arbitrary code 
source, and let M_ be the free monoid generated by this alphabet. 
The sequences emitted by this source are, then, words in the free 
monoid M,. To say that the given source is ergodic is, then, to 
say that for each sufficiently long word w emitted by the source, 
we must have 


A, (w) 
A(w) 


= Pi t= 1,eeeyn, 


where \,, A are the length mappings defined on p. 74, and p, isa 
constant. From this definition, it follows immediately that 


n 


Y= 1 


1 


Conversely, let there be given n numbers p,, P,,+++,p,, Satisfy- 


n 
ing the condition nS p, = 1 (for convenience we restrict ourselves 


im1 
to rational numbers p,).' Consider the totality of all words weM, 
such that for each 2, 


d, (w) 
d(u) a Pi 
and denote this set of words by Eo 3 We prove 


Theorem : E. a is a submonoid of M_. 

Proof. Let w,, wy be words in Eu > we must show that w,w, 
is likewise a word in E,, Since the mappings A;, A are miene 
to be monoid homomorphisms of M, into Z*, we have, because the 
composition on Z* is ordinary pues 


A, (uw, v2) z A, (w,) +A, (w,) 2) 
A(w, Wg) e A(w,) + A(uw,) 


86 ROBERT ROSEN 


But by hypothesis, A,(w,)=p,A(w,) and ,(w.) = P; A(w,) for 
taal yoo 6 Me Substituting these expressions in (2), we find that 


A ,(w, wv.) 


= , l= 1 © oo @ the 
: : i? 3 ’ 
A(w, w,) 


Thus, w, “,€ EY i If we now adopt the convention that the empty 
word e dabetion ‘the condition A, (e)/A(e) =p, for each z, ‘then it 
follows—that ter ) is a submonoid iat Ms as Peeertea: 

Hence, it is seen that to specify an “ergodic submonoid of a free 
monoid M, is equivalent to specifying an ergodic code source 
using the set of generators of M, as an alphabet. . The importance 
of ergodic coding lies in the fact that virtually all natural and 
man-made coding procedures are of this type; for example, all 
modern languages, such as English, can be characterized as being 
of an ergodic nature in the above sense [for greater details on this 
matter, see, C. Shannon and W. Weaver (1949)]. Hence, it is 
natural to investigate whether the DNA-protein coding is of an 
ergodic character.. To simplify matters, we shall assume at the 
outset that actual DNA molecules are sufficiently long so that 
any Statistical properties of the code will be found in each mole- 
cule, with a probability which is negligibly different from unity. 
Since experiment indicates that the order of magnitude of the 
length of natural DNA molecules is of the order of 2 x 10° (Beadle, 
1957), this assumption is probably sufficiently accurate. 

Let us now return to the main line of our development... We 
shall term the submonoid E - which we have constructed above 
the ergodic submonoid of M, determined by the given numbers 
Pyrte+)P,+' It follows from Theorem 3 that any ergodic sub- 
monoid of a free monoid is itself free. However, we now prove a 
rather interesting result: 

Theorem 5: An ergodic submonoid of a finitely generated free 
monoid is free on a countably infinite set of generators. 

Proof. Let the symbols {a,, @,,+++,@,} generate the free 
monoid M_, and let the numbers {p,, Por*++,P,,$ determine the 
ergodic submonoid E pi)? where 0 <p, <1 for each 7%. We may 
put the numbers p; over ‘a common denominator by writing p,; = 1,/r, 
where ris chosen as small as possible. We now observe that all 
ergodic words of M_, i-e., words of M_ which lie in E. p;) are 
necessarily of length nr, where n is some non-negative sf Se 
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Let us agree to call an ergodic word prime if it cannot be ex- 
pressed as the product of ergodic words of smaller length. Every 
prime word is necessarily a generator of F p;) hence, our theo- 
rem will be proved if we can show that for every integer n, we 
can find a prime word of length nr. 

Let an arbitrary ergodic word of length nr be divided into seg- 
ments of length r, as below: 


iste. «ca tl. | @, "tk ee ae ee : wy 
ty tg tr tr+y try tor inr—r lUnr 


1 9 nn. 


We shall show how to construct a prime word from the given word. 
Let us consider first segment 1 of the displayed word. If this 
segment is not ergodic, we leave it alone. If it is ergodic, then 
we may. remove the ergodicity as follows: we replace any one of 
the symbols which occur in segment 1 with the first different 
symbol which we encounter in the terminal segment of length 
(r-—1)n of the displayed word. This procedure does not affect 
the ergodicity of the whole word, but gives us a new ergodic 
word w’, in which the initial segment of length r is not ergodic. 

Next we consider the segment 1+ 2 of length 2r of the word 
u’. If this segment is not ergodic, we once again leave it alone. 
If it is ergodic, we can remove the ergodicity by removing any 
symbol of segment 2 and replacing it by the first different symbol 
in the terminal segment of length (r— 2)n of the word uw’ Once 
again it is clear that by means of this permutation of symbols we 
obtain a new ergodic word w’’ whose initial segments of length r 
and 27 are not ergodic. 

By proceeding in this manner, we are able to construct an 
ergodic word with the property that none of its first (r—1)n seg- 
ments are ergodic. Hence, this word is prime, and the theorem is 
proved. 

As a concrete illustration of the technique used in the proof of 
Theorem 5, let us choose m = 2, p, = 4, Pox =. ‘Here r= 3; itis 
evident that E word w of M, cannot satisfy the conditions 
rd (w)/A(w) = = A 9 (u)/r (w) = = 2 unless w has length 3n for some 
sates Ne . The re of length °s in the ergodic submonoid of M, 
determined by (4; =) are 


A1%q_% G0, Vs ober hee 
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We shall construct an ergodic word of length 2-3 which is not a 
product of these words. . Consider the word 


UW=4,4,% | Ay Gy, 


which is an ergodic word of length 6 expressible as a product of 
ergodic words of length 3.. To obtain a word which is not thus 
expressible, we follow the procedure outlined in the proof of 
Theorem 5.. Explicitly, take the symbol a,, ‘which appears in the 
first segment of w, ‘and interchange it with the first symbol a, of 
the second segment of uw. We thus obtain a word w’, 


which is still ergodic, but which is obviously not a product of 
ergodic words of smaller length; i.e., w’ is prime. 

In order to apply Theorem 5 to our coding problem, we must first 
investigate the notion of ergodicity we have defined in the case 
when the free monoid containing the ergodic submonoid is itself 
embedded into a third free monoid. More explicitly, let us sup- 
pose that M_ is generated by the symbols {b,,'},,..., 6,3, and 
that N is an ergodic submonoid of M,. By definition, this means 
that we are given m rational numbers p, with sum unity such that 
A, (w)/A(w) = p;, weN. By virtue of Theorem 2, we know that 
M,, can be embedded as a submonoid of the free monoid Mo 
generated by a pair of symbols {a,, a,}.. The generators of M., 
when embedded thus will obviously be a set of words {uw 
Wo, ster, Uw } Of My. 

We now encounter the following problem: M, has its own length 
mapping, which we shall denote by yp, ‘and which is defined as 
follows: p(w) = the number of symbols which occur in the word 
weM,. Likewise, M, has its own.mappings p,,‘p,, which are 
defined as the number of times that a, and a,, respectively, ap- 
pear in a word of M,-' Now N was a submonoid of M_, ergodic 
with respect to the mappings A,, A, (¢=1,...m) which were in- 
digenous to M_.' Under the given embedding of M_, the sub- 
monoid N now becomes a submonoid of M,. We may, therefore, 
inquire whether the submonoid N remains. ergodic (with respect 
to new rational numbers 7,, 7,, Of course) when referred to the 
mappings p, #,, #4. The answer is in the affirmative: 


1? 
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Theorem 6: If a submonoid N CM_, is ergodic, it remains ergodic 


under an embedding of M_ into M,,. | 
Proof. As above, let us suppose that M_ is generated in M, by 


the words {w,, 


Wo +++, w,_}.' Let us agree to write p,(w,) = 7, 


p,(w,) = s,-' Then, for any word weM,, we have 


3 ri A, (w) = p(w), 
zs (3) 
> 8:4, (x) = pL, (wv). 


2. 


If N is an ergodic submonoid of M_, we have for all words weN, 


A (w) = p,A(w); 4 m-Lyceey Me 


Substituting this expression in (3), we find 


Therefore, we 


1 (0) = Jo rp, Mw) = we) 7 10,5 


owe | i= 1 
pea(w) =) 8,p,d(w) = A(w) D1 8,2; 
t= P= 


have, since p(w) = p,(w) + p,(w), 


d(w) >> tP; 
pe, (w) int 


(w) + p(w) m n 
p(w) + pale aw) (Yo raeit Do sar 


and likewise, 


i= 1 i=1 
xr.p. 
X(r, + 8;)D; 
[Ly (w) x 85D; = ( 


———_____—_—_———_ = = TT 5 
p(w) + pg(w) Er, +8,)r,  * 
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It is evident that 7, + 7, = 1, ‘and hence, ‘by definition, N is an 
ergodic submonoid of M,. This proves the Theorem. . 

Theorem 6 is the simplest possible situation. However, exactly 
the same line of argument yields 

Theorem 6a: An ergodic submonoid N of a free monoid M_ re- 
mains ergodic (in the sense above) in any free monoid M, in which 
M_, 1s embedded. 

Proof. Left to the reader. 

We are now able to draw some interesting biological conclusions, 
based on Theorems 5 and 6 (or better, 6a). First, let us observe 
that the general DNA-protein coding problem, as formulated above, 
reduced to the problem of embedding the free monoid M,, zso- 
morphically in the free monoid M,. However, it follows from 
Theorem 5 that M,,, when so embedded, is no¢ an ergodic sub- 
monoid of M,, ‘since M,, is, by definition, free on only a finite 
set of generators.. There seem to be only two possible ways of 
explaining this circumstance: 

1. The DNA-protein code is not ergodic. 

2. The DNA-protein code is ergodic, but not every conceivable 
polypeptide can be coded, though all may exist in nature; those 
which can must, then, form an ergodic submonoid of Moo? and 
hence, of M, (by Theorem 6a). 

If the first alternative holds, then we are faced with the un- 
likely situation that one of the fundamental coding processes in 
nature deviates from the simple statistical structure which one 
would tend strongly to assume that it possessed.. The second 
alternative seems more satisfactory in this regard, and although 
it is more restrictive than one might like, it may perhaps provide 
a theoretical justification for the fact that only a small fraction of 
the possible polypeptides have been observed in nature. . However, 
it seems very difficult to decide experimentally which of these 
alternatives actually holds. 

The seeming fact that one of the two possibilities noted above 
must be satisfied gives us interesting and not altogether obvious 
suggestions about the properties of the DNA-protein code. We 
may now go one step further, and ask if it is possible for our 
analysis to throw some light on the actual structure of the code 
itself. We shall now attempt to show that the machinery of Theo- 
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rem 6 enables us to give at least a partial solution of the struc- 
ture of this code. 

Let us first set up this machinery for the case of the special 
embedding (of M,, into M,) which we have shown to be the alge- 
braic formulation of the coding problem. We are given the free 
monoid M,, ‘generated by the set {a, ,..,@,}, and the mappings 
By p,tM,— Z* (i =1,...,4), where (2) denotes the number of 


4 
times the generator a, appears in w, and p(w) = DH; ()- We 
i=1 

then assume that M,, is embedded in M, and generated by the set 

of words LWyyeees Wools w,eM,; we further assume that the as- 

sociated mappings M,, > 7* are denoted by A, A, (j= 1,+-- , 20). 

If now N is an ergodic submonoid of M,,,, relative to the rational 
20 

numbers {p,,+++)Pg.} (where 0< Pp, < 1 and >; = 1), then the 
j= 

technique of proof of Theorem 6 gives that N is an ergodic sub- 

monoid of M,, determined by the rational numbers bajyeee'gmyhs 

where 


Sy 00105 ee oe en eee eet a ee (4) 


4 
py Se ye oP; 
tt 


i=1 j=1 


and we have written ii =p, (w,). 

We now inquire into the biological significance of the symbols 
appearing in the set of equations (4).' The numbers 7, ‘p, (i= 
1,%..,4; j=1,-...,20) define the statistical structure of the 
coding process by giving the relative frequency of occurrence of 
the generators of M, (which we may call nucleotides) and M,, 
(amino acids), respectively, in the words of the submonoid N (the 
codable polypeptides). Stated otherwise, if w is a codable poly- 
peptide, ‘then the 7, determine the composition (in terms of relative 
amounts of nucleotides) of the DNA molecule which codes w, 
while the P; determine the composition (in terms of relative 


92 ROBERT ROSEN 


amounts of amino acids) of the polypeptide w itself. In particular, 
for each fixed j = 7,, ‘the numbers Be (i =1,...,4) determine the 
relative amounts of nucleotides occurring in the word Ul 5 which is 
one of the generators of M,, under the given code; i.e., these 
numbers give the percentage composition of the nucleotide chain 
which corresponds to the individual amino acid w, i bi 

Now the numbers 7,, p, can, in principle, be determined empiri- 
cally, by investigating the percentage composition of DNA mole- 
cules and polypeptides; these should be constant for all su/f- 
ficiently long molecules. We may then employ the set of equations 
(4), which consists of four equations in the 80 unknowns 7‘. We 
can solve these equations for any four of the eee thereby obtain- 
ing at least some concrete information about the actual code 
structure (subject to our hypotheses). 

A fuller determination of the structure of the DNA-protein code 
is, then, seen to be dependent upon obtaining more relations be- 
tween the 7‘. Our search for further relations has thus far been 
unsuccessful, but continued investigation in this area should 
throw more light upon the details of the code structure. 

Since we have not succeeded above in determining the DNA- 
protein code uniquely, the question naturally arises as to whether 
there exists a unique solution to the problem; i.e.,:we may ask 
whether there can exist more than one DNA-protein code in nature 
satisfying the condition of ergodicity. Although we cannot at the 
moment supply a full answer to this question, we shall supply a 
brief discussion of certain relevant aspects of this problem. 

Let us first consider a problem which arose in pure mathematics 
around 1900.' Let @,,'a,,+.., a, be symbols which generate the 
free monoid M_.' A word veM, is said to be embedded in a word 
eM if there exist words Wy, ‘Wy, in M such that #= U, Wey 
Now let there be given a finite subset of the cartesian product 
M,, xM_3 i.e., a set of pairs (U1, %1)) (Ug, Ug), +++, (U,, ¥,), Where 
wu, v,eM_ ,2=1, 2,...,n. Such a subset is sometimes referred 
to as a dictionary. Two words a2, yeM, are called immediately 
equivalent under a given dictionary if they are, respectively, of 
the form 


T= U,U;W,, Y=W,%;,U, 


LENS 
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for some i, 1 Si <n, and some words w,, w,eM,. Two words 
x, y are called equivalent if we can find a finite sequence of 
words a, #,,-++,%, such thatz, =a, 7 =y, and a, is immediately 
equivalent to z,_, for? =2,...,7 

The problem mentioned above consists of finding an algorithm 
(i.e., a decision procedure terminating in a finite number of steps) 
for deciding whether two given words are equivalent under a given 
dictionary; this question is known as Thue’s Word Problem. This 
problem, however, is known to be unsolvable in a certain well- 
defined sense; for full details see, e.g., Kleene (1952). 

To return to the main argument, let us now suppose that there 
are two distinct DNA-protein codes in nature. According to our 
initial abstract formulation, this would mean that there are two 
distinct monomorphisms f, g:M,, — M,- IfM,, is generated by 
the symbols {a,, @ 


re then we can construct a dictionary 


2? ao!s 
in M, by considering the set of pairs {[f(q,), g(a,)\, [fCa,)s 
g(a,)I, oeesIf(Gy5)s g(a) J}. The physical meaning of this dic- 
tionary is that two DNA molecules will be considered equivalent 
if, and only if, they code the same polypeptides in the respective 
codes. But the unsolvability of Thue’s problem means that there 
may be no finite logical procedure for determining such an equiva- 
lence, although there is a definite physical procedure for doing 
this (namely, comparing the polypeptides which result from the 
biological action of the two DNA molecules which are to be 
tested for equivalence). It may be considered a contradiction to 
have a decision process which may be carried out physically, but 
yet is not describable by means of a finite mathematical proce- 
dure. Hence, the conclusion is strongly suggested that if every 
polypeptide is capable of being coded by DNA, there cannot exist 
two distinct DNA-protein codes. 

If the actual DNA-protein coding processes are ergodic, how- 
ever, then not every polypeptide can be coded, and the question 
as to whether more than one such code exists in nature is not 
settled by the above argument. In this case, we have two mono- 
morphisms A:N —> M,, k:N’ + M,, where N and N’ are ergodic 
submonoids of M,,. However, we must have N and N’ disjoint in 
Pe otherwise, the same argument as we have just employed 
would set up a dictionary in M, as before. 
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If this disjointness condition is satisfied, then, there seems to 
be no a priori reason why more than one ergodic DNA-protein code 
cannot exist in nature. In fact, the existence of more than one 
code may be indicated by the experimental determination that the 
ratios of the generators of M, (i.e., the polynucleotide ratios) vary 
from species to species (Chargaff, 1955). This may be explained 
either by the presence of redundant polynucleotides, or by assign- 
ing to each species its own ergodic code. Under this last hy- 
pothesis, the disjointness condition implies precisely that no 
two species can have a protein in common. That this actually 
seems to be the rule in biological systems is indicated by nu- 
merous antigen-antibody experiments, and in a more direct fashion, 
by the direct determination of amino acid sequences in the homolo- 
gous proteins of various species (in particular, Sanger’s work on 
the insulin molecule).' However, these last remarks must remain 
in the nature of speculation until further experimental work is 
forthcoming. 


This work was aided by United States Public Health Service 
Grant RG-5181. 
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A NOTE ON TOPOLOGICAL BIOLOGY 
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In line with a recent suggestion by the author (Bull. Math. Biophysics, 
20, 267-73, September, 1958) that not only does the organism as a whole 
map on the primordial, but that each organ can also be thus mapped, it is 
shown that the previously introduced abstract spaces, which represent an 
organism, contain subspaces which map continuously on the space of the 
primordial. Several theorems about those subspaces are proven. 


In a previous paper (Rashevsky, 1958b) we suggested a possible 
extension of the principle of biotopological mapping (Rashevsky 
1954; 1958a; 1959) by postulating that not only the organism as a 
whole could be mapped on the primordial, but that each organ could 
also be mapped on the primordial organism. This leads to the im- 
mediate conclusion that the whole organism can be mapped on each 
of its organs. 

In terms of the set-topological representation of the organism, 
discussed previously (Rashevsky, 1958a), the above postulate 
amounts to stating that in the abstract space S,, which represents 
the organism, there exist subspaces S,*, such that S, can be con- 
tinuously mapped on each of the S,* and that every point of S5* 
maps on itself. Such subspaces, if they are compact subspaces in 
a Euclidian space FR”, are called retracts (P. Alexandroff and 
H. Hopf, 1935). Inasmuch as our abstract space S, is anything but 
Euclidian, the terminology cannot be used here, lest it cause con- 
fusion. We shall conveniently call the subspaces S* of S, which 
have the above property, generalized retracts. ' 

As we have seen (Rashevsky, 1958b), the requirement that the 


organism maps on each of its organs leads to some conclusions 
which are verifiable by experiment. The discussion of the spaces 


So (Rashevsky, 1958a) also has led to a number of conclusions 
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which agree with known facts. It is, therefore, important to inves- 
tigate whether the above extension of the principle of biotopologi- 
cal mapping imposes any restrictions on S,, or even is incompati- 
ble with the properties of S,.° We shall see that actually any Se 
does contain a large tens of subspaces S,* and that, cheretorel 
the suggested extension is quite compatible with the results ob- 
tained previously. 

We shall use the same notations as before (Rashevsky, 1958a) 
and shall first prove the following: 

Theorem 1. Let n, be the number of different subproperties P, , 
in S, which are included in P, in S,. Then the number of (not 
necessarily disjoined) generalized retracts in S, is greater than 
the largest n,. 

Proof: Let the largest n, correspond to the particular z = 7*, 
Take an arbitrary Ping amongst those n, subsets. If P, —> P, in 


Sp, then we shall ete Pit og al gape In other words, in the graph 


associated with S,, the point Px. og will be connected by an arrow 
to the point Pi og? Tf, in its turn, P, — P, in Sp, then Pio, meytt ST 


and so forth. Thus, starting with P. we construct a directed 


*X 
graph Cin cg , which is homeomorph to the directed graph of S,. 
Since (Rashevsky, 1958a) the directed edges of the graph define the 
neighborhoods of the points in S,, therefore, the vertices of Gro, 
form a subspace of S,, which is peeataniac to Sp. Since S, maps 
continuously on S,, therefore, it also maps continuously a the 
subspace which we just constructed, and each point of this sub- 
space maps onto itself. Hence, this subspace is a generalized re- 
tract of S,. 

We may use the same procedure starting with any of the n « points 
in S,, and each time obtain a different generalized retract of Sige 
Thus, the number of those generalized retracts is at least as la 
ae 

To prove that it is actually greater than n,,, we remark that, be- 
cause of the connectedness of S,, there exist such a & and such a, 
and 4, that P,, — P,, and Pr on — Prop where 4 is different 
from 4,, if P, P, in s, eae: 1958a, p. 80), 

But Pro, _ Page if ’P, — P| in Sp. Similarly, Pra, — Proy* 
In general, 4, is dittogont from 4. Starting again with Px» fol- 
lowing until the point i o, but passing from there to Pig, instead of 


Pre. ,we again obtain a a Civ cr? which is homeomorph with the 
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graph of S,, but is different, at least in some points, fromG,, . 
This graph Gis cg again defines a generalized retract of S,, which 
is different from the one obtained previously.’ Hence, to each of 
the n, points there corresponds more than one generalized retract. 
This proves the theorem. 

Theorem 2. If S, is maximal (Rashevsky, 1958a, p. 79), then the 
number of different generalized retracts in it is equal to 1g n,, taken 

t 


over all the possible indices 2. 

Proof. Choose any point P,,. If P; — P, in Sp, then the total 
number of arrows which connect any of the n, points P,, to any of 
the n, points P, ,is njn,. If P, — P, in Sp and if S, is maximal, 
then every one of the n, points is connected to every one of then, 
points. In the graph which is associated with S,, the total number 
of directed ways from any of the n, points to any possible point of 
the n, points is 7,7); and since there are n, points P,., therefore, 
the total number of directed ways from all of the n; points to all of 
the n, points is nn,n, If P, — P, and P, — P_, in Sp, then for 
each way between any of the n, points P,,, and any of the 7, points 
P,,, we can add another n,, branches, thus obtaining a total of 
nn, nn, different, though not disjoined, trees of the form: 


The same is readily seen to hold when we have a uniformly or 
non-uniformly directed cycle (Rashevsky, 1955). Hence, continuing 
the procedure until all indices ¢ are used up, we find the the total 
number of possible graphs IT n,. Therefore, this is also thenum- 


ber of the generalized retracts. 

Theorem 3. The maximum possible number of disjoined general- 
ized retracts in S, is equal to the smallest n,. 

Proof. et the smallest n, correspond to 7=72’. Since every 
generalized retract must map continuously on Sp, therefore, it must 
contain a point which maps on P,, in other words, one of the points 
Pv ye If there were more than n, disjoined generalized retracts, 
some of them would have at least one point P,,, in common and, 
therefore, would not be disjoined. 

Whether or not the above theorems have any biological applica- 
tions cannot be said at present.’ However, in line with the general 
ideas suggested elsewhere (Rashevsky, 1955), a study in abstracto 
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of the properties of spaces which represent organisms must of ne- 
cessity precede any applications. By studying the mathematical 
properties of such spaces like S,, we may hope to find some prop- 
erties which, when translated into biological terminology, will de- 
scribe known biological characteristics of the organism. 

Very tentatively, we may suggest the following: 

It does not seem to be biologically meaningful to consider every 
generalized retract of S, as describing the biological relations in 
organs of the organism represented by S,. It may, however, be 
possible to assume that the total number of organs in a fully dif- 
ferentiated organism is equal to the total number of disjoined gen- 
eralized retracts.’ In biological terms this will mean that the maxi- 
mum number of distinct organs that an organism can have is equal 
to the smallest number of differentiated subproperties of any of the 
basic biological properties of the organism. It may, perhaps, be 
more plausible to identify the ‘‘generalized retracts,’’ not with 
organs, but with cell types. The interesting problem then arises of 
finding such sets of ‘‘generalized retracts,’? which correspond to 
organs. This will be discussed later, » 


The author is indebted to Mr. Robert Rosen for a discussion of 
this note, 

This work was aided by United States Public Health Service 
Grant RG-5181, 
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In a preceding paper (Bull. Math. Biophysics 20, 71-93, 1958) the prin- 
ciple of biotopological mapping was formulated in terms of a continuous 
mapping of an abstract space, made from the set of biological properties 
which characterize the organism, by an appropriate definition of neighbor- 
hoods. In this paper it is shown that we may consider directly the map- 
pings of the different sets of properties which characterize different 
organisms without taking recourse to abstract spaces. All the verifi- 
able conclusions made in the preceding paper remain valid. A serious 
difficulty mentioned previously is, however, avoided and the possibility 
of more general predictions is established. 


In a series of preceding papers (Rashevsky, 1954; 1955a; 1955b; 
1955c; 1955d; 1956a; 1956b; 1956c; 1957a; 1957b), a new biologi- 
cal principle, called the principle of biotopological mapping was 
proposed and developed on the basis of the theory of graphs. In 
another paper (Rashevsky 1958a), hereinafter referred to as loc. 
cit., a set-topological approach was introduced and shown to have 
a number of advantages over the previous approach. Subsequently 
(Rashevsky 1958b), a connection between the two approaches was 
established. For the understanding of the following, familiarity 
with loc. cit. is presupposed. The notations used here are the 
same as in loc. cit. 

In loc. cit. we constructed a topological space which corresponds 


‘to an organism, by defining the neighborhoods in the set of biologi- 


cal properties in terms of the relation of immediate precedence or 
immediate succession of two biological properties. On page 91 of 
loc. cit. we pointed out that a difficulty of this approach, a di ffi- 
culty which is also ch aracteristic of the graph-theoretical approach, 
consists in the circumstance that two biological properties, of 
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which one immediately precedes the other, are usually connected 
also by a number of other relations which may he quite different 
from the relation of immediate precedence. 

We may, however, avoid the construction of a topological space 
and speak only of sets of hiological properties. As we have seen 
in loc. cit., those properties of an organism which are subproper- 
ties (or subsets) of a given property of the primordial, map on that 
given property in the primordial. If we wish the mapping to pre- 
serve the hasic relation of immediate precedence and immediate 
succession, we do not need to introduce the spaces S, and Sp 
(loc. cit.), but may operate directly with the corresponding sets 
Si and Sf (loc. cit., p. 77). All we need is to make the following 
requirement: 

(A’). If in Sf the biological property P, immediately precedes 
the biological property P,, then in Sf at least one subproperty 
pee immediately precedes at least one suhproperty 69 me If in 
Sp the biological property P; does not immediately precede P,, 
then no P,, immediately precedes any P, 2. 

This requirement takes the place of the postulate (A) in loc. citt. 

In loc. cit., however, we also introduced the postulate (8) which 
requires that the space S, be connected. As we have seen, this re- 
quirement is necessary, hecause without it the possibility would 
exist that a multicellular organism is characterized by a collection 
of sets of biological properties, such that there is no interaction 
whatsoever between those sets, This would be a biological ab- 
surdity. Thus, the connectedness of the abstract space S, reflects 
not so much the geometrical connectedness of the organism as its 
‘‘physiological connectedness.’’ Geometrically, strictly speaking, 
a multicellular organism is not connected, The erithrocytes and 
leucocytes in the blood of a metazoan represent geometrically non- 
connected living parts of the organism inasmuch as the plasma in 
which they are suspended is itself considered to be non-living. 
But, physiologically, those cells are definitely connected to other 
parts of the organism. ° 

Now that we speak only in terms of the set S%, rather than of the 
space S|, we cannot as conveniently formulate the requirement of 
physiological connectedness, as we did in the postulate (B), (loc. 
cit., p. 79). Herein lies an advantage of describing the organism in 
terms of a space, rather than in terms of a set. We now must sub- 
stitute for (B) the following more lengthy requirement: 
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(B’). If in S;, the biological property P; immediately precedes P,, 
then for at least one 7 there exists a P.,, in 8S,” which immediately 
precedes at least two Pig's in S’. 

That this requirement assures a physiological connectedness of 
the organism is seen by an argument which is patterned after the 
argument on page 80 of loc. cit. We may again consider n, =n, 
(loc. cit. p. 79). Then in absence of requirement (8’) we would ob- 
tain n sets of subsets P., (i and « variable), such that no subset 
of one set immediately precedes a subset of another set. The or 
ganism would be physiologically non-connected. 

Postulates (C) and (D) of loc. cit. (page 81) remain unchanged 
except for the substitution of the terms “‘sets S*,’? “‘set 8,” 
and ‘‘S’*? for ‘‘spaces S,,,”” ‘‘space S,,’’ and “‘S,,’’ correspond- 
ingly. 

It is seen that all conclusions drawn in loc. cit. remain valid 
with this new formulation. But the latter gives us more than that. 
Instead of the relation of immediate succession and its converse, 
immediate precedence, we may now consider any kind of binary re- 
lation between P. and P,. We shall designate it by 


Pi Ry, Pye (1) 


We write R,, in order to indicate that while between a given P, 
and a given P, one particular relation may hold, for another pair of 
P, and P,, in other words for other values of z and &, the relation 
may be quite different. In general, several different binary rela- 
tions will hold between a given P, and P,. Thus, F,, is, in gen- 
eral, the intersection of several relations R, R®,..., i>. 


For example, it may happen that a property P, not only immedi- 
ately precedes the property P,, but also a sudden release of a pro- 
longed inhibition of P, results in an increase of the intensity of 
the property P,.° If we denote the relation of immediate precedence 
by R®, and the just mentioned, other relation by R®, then we have 
PRP, with RB, = RO o R®,: On the other hand, another property 
P_ may not only precede the property P,, but a short, complete in- 
hibition of P, results in a long inhibition of P,. Let us denote this 
relation by R“. Moreover, as an example, it may happen than an 
increase in the intensity of P, does not result in an increase of in- 
tensity of P,. Denote this relation by R‘®, We now have P_R,.P,, 
with Rp. = RPADROOR®, 
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Now instead of (A’) and (B’) we introduce the following postu- 
lates: 

(A,). If and only if in Sf we have P;F,, P,, then in S* there ex- 
ist at least one P,, and one P, 2 such that P;, h;,P, 2 and for at 
least one 7 we have P,,R;, Py g and P;,R;,P,. (B # y)s 

Postulate (A,) replaces (A’) and (B’). It assures that if all 
P,,7s map on P,, then all relations are preserved under such a 
mapping. 

By making F,, = R® = relation of immediate succession, we re- 
duce the situation to that discussed previously. Hence, again, all 
conclusions drawn in loc. cit. remain valid,’ But to the extent that 
R., now stands for any binary relation, we may expect a number of 
other more general and interesting conclusions to be drawn. This 
will be investigated elsewhere. The method of derivation is basi- 
cally the same as in Joc. cit. If in an organism we haveP;,, F;, P, a, 
then it follows from (A,) that in the primordial we have P.R,, P,. 
Hence, from (D) (loc. cit.), it follows that there exist organisms, in 
which P,F,,P,,. From (C) (loc. cit.) it then follows that there 
exist unicellular organisms in which P;,f;, P;,, unless the ex- 
istence of either P, , or P;,, or of the relation F,,, is incompatible 
with the size of a single cell. 

In biological systems we are, however, not limited to binary re- 
lations. There may be m-ary relations between any m properties. 
As an example of a possible ternary relation (m = 3) we may men- 
tion the following: The biological property P, is manifested if and 
only if both biological properties P, and P, are of sufficient inten- 
sity.’ Let us denote this relation by R;,,(P,,P.,P,). The study of 
the preservation of such type of m-ary relations under the biologi- 
cal mapping is likely to lead to much more important results than 
the previous case of binary relations. Denoting by 2 awe Ba 6 a : 
P;,++.P;) an m-ary relation, we now require ; 

A,). If and only if in Sf we have Rij,...im Pe PE then in 
S* there exist ¢wo or more sets Pia, Pr Py Om? for the same set 
of indices ¢,,..., ¢,,, but for different sets of Ay, eee, X, Such 
that R; reeeg (Pi og? Pere P, tm) holds. 

The difficulty inherent in the graph-theoretical approach and in 
loc. cit., which we have now overcome, has already been overcome 
in an ingenious manner by Robert Rosen (1958, 1959), who intro- 
duced the notion of an ‘‘abstract block diagram,’’ or ‘‘inverted 
graph,’’ using the theory of natural equivalences. It is now neces- 


“s 


sO 


oil atk) So 


PIPE eh 


x 


SET-THEORETICAL APPROACH TO BIOLOGY 105 


sary to investigate the possible connection of Rosen’s approach 
and ours. ° 

The requirement ‘‘two or more’’ in (A,) insures again the ‘‘phys- 
iological connectedness’’ of the organism. However, while the 
corresponding formulation in (A,) is the simplest conceivable, 
there are now several possibilities of equal simplicity. The one 
given in (A,) is, therefore, to be considered only as tentative. It 
has been chosen because it leads directly to an argument similar 
to the one used previously.’ Namely if in an he for a given 
set 7,, «2+, %,, we have R;, .s.;) (Piay---s P; aes ), then it 
follows from (A,) that in the farimardial we have RF; , ...; (P. 
P; de eee there exist organisms in which we > have R; 


eae - Pee) where p, # & (p= 1,2,..., m)e 
With this generalization of our previous results, it is no longer 
appropriate to speak of the principle of biotopological mapping. 


The designation ‘‘principle of biological mapping”’ is suggested, 
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The author is indebted to Mr. Robert Rosen for a discussion of 


this note. 
This work was aided by United States Public Health Service 
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